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SECTION  1 


OVERVIEW,  CONCLUSIONS  AND  RECOMMENDATIONS 


OVERVIEW 


Deep  based  facilities  will  almost  surely  be  deployed  well 
beneath  the  water  table.  In  addition  to  the  usual  problems  with 
inflow  of  water  during  construction,  such  facilities  may  be 
subjected  to  additional  water  problems  resulting  from  the  high 
intensity  dynamic  loadings  expected  during  a  nuclear  attack.  At 
the  Nevada  Test  Site  (NTS)  test  structures  adjacent  to  contained 
nuclear  bursts  have  sustained  only  minimal  structural  damage  from 
the  dynamic  stresses  and  ground  motions,  but  have  been  completely 
flooded  with  ground  water  mobilized  by  the  blast.  Similar 
occurrences  at  key  locations  on  a  deep  based  strategic  system 
would  defeat  the  system  even  though  it  survived  the  dynamic 
threat. 


This  study  identifies  the  principal  mechanisms  respon¬ 
sible  for  explosively  generated  inflow  of  water  and  makes 
preliminary  estimates  of  the  potential  magnitude  of  this  problem 
at  a  typical  deep  basing  site  as  exemplified  by  Generic  Mountain  C. 
A  summary  of  hydrologic  response  data  from  NTS  underground 
nuclear  and  high  explosive  (HE)  tests  beneath  Rainier  Mesa  is 
used  to  identify  similar  response  mechanisms  which  would  threaten 
deep  based  facilities  subjected  to  surface  or  shallow  buried 
nuclear  bursts.  Two  phase  one  dimensional  axisymmetric  code 
calculations  are  used  to  assess  the  potential  threat  from 
explosively  generated  inflow  of  water  at  Generic  Mountain  C. 

Flow  rates,  flow  volumes  and  pore  pressure  dissipation  for  a  wide 
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variety  of  rock  types  ranging  from  the  impermeable  tuffs  of 
Rainier  Mesa  to  the  permeable  sandstones  of  Generic  Mountain  C 
are  calculated  using  the  CONSL  two  phase  code. 

The  NTS  data  base  is  summarized  and  highlighted  in 
Sections  2,  3 ‘and  4.  Section  2  describes  the  unique  geologic  and 
hydrologic  regime  of  Rainier  Mesa,  with  its  high  porosity,  very 
low  permeability  altered  tuff  and  variable  perched  water  table. 
Section  3  summarizes  the  most  significant  and  best  documented 
instances  of  inflows  encountered  during  construction  of,  and 
exploration  for,  the  E,  N  and  T  tunnel  complexes  beneath  Rainier 
Mesa.  Inflow  data  from  exploration  and  construction  activities 
are  instructive  in  identifying  potential  water  sources  and  in 
describing  interconnectivity  and  extent  of  water  bearing 
fractures.  The  fourth  section  presents  data  illustrating  the  two 
explosive  induced  water  producing  mechanisms  identified  at  NTS. 

The  subsection  on  the  MIGHTY  EPIC  event  describes  the  mobil¬ 
ization  of  fracture  water  by  a  contained  nuclear  detonation  and 
documents  the  flooding  of  test  structures  and  access  drifts  which 
occurred  on  that  event.  The  subsection  on  the  ONETON  HE 
contained  burst  describes  the  production  of  water  from  the  residual 
stress  field  produced  by  cavity  expansion.  The  residual  stresses 
squeezed  water  from  the  pores  of  the  rock  over  a  period  of  many 
days.  Migration  of  the  pore  water  away  from  the  higher  stress 
regions  near  the  explosively  formed  cavity  was  accompanied  by 
pore  collapse  and  a  slow  decrease  in  residual  stress. 

Sections  5  through  8  investigate  the  mobilization  of  pore 
water  by  residual  stresses  and  pore  pressures  and  use  parametric 
two  phase  flow  calculations  to  estimate  the  magnitude  of 
potential  flow  due  to  this  flow  mechanism  for  a  deep  based  system 
beneath  Generic  Mountain  C.  Section  5  and  Appendices  A  through  D 
describe  the  CONSL  code  used  in  the  parametric  calculations. 


3 


f 


Section  6  describes  mechanisms  by  which  residual  stresses  and 
pore  pressures  are  developed  from  both  contained  and  surface 
bursts  in  or  on  saturated  rock.  The  parametric  flow  calculations 
are  presented  in  Section  7  for  a  broad  range  of  rock  types 
encompassing  the  NTS  tuff  and  the  rocks  making  up  Generic 
Mountain  C.  Sensitivity  studies  of  flow  volume  as  a  function  of 
changes  in  rock  properties  are  also  presented.  Finally,  a 
calculation  of  flow  on  the  ONETON  event  using  measured  residual 
stresses  and  rock  properties  is  presented  in  Section  8  and 
compared  to  the  ONETON  data.  The  agreement  between  the 
calculation  and  the  data  validates  the  calculat ional  procedures 
and  strengthens  the  concerns  raised  by  the  calculations  of 
Generic  Mountain  C  in  Section  7. 


CONCLUSIONS 

•  Post  shot  water  has  been  produced  by  two  mechanisms  at 
NTS. 

Complete  flooding  of  test  structures  occurred  on 
the  MIGHTY  EPIC  nuclear  event.  Flooding  was  caused 
by  fracture  water  mobilized  by  the  dynamic  stresses 
and  ground  motions  (including  block  motion).  Water 
entered  the  structures  through  relatively  minor 
tears  in  the  steel  liners. 

—  Flooding  of  the  instrumentation  drift  on  the  ONETON 
HE  event  resulted  from  interstitial  pore  water  being 
driven  out  of  the  voids  by  the  residual  stress 
buildup  around  the  cavity  formed  by  the  explosion. 
Such  flooding  by  pore  water  migration  has  not  been 
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a  major  problem  in  the  NTS  tuff  because  of  the  tutf's 
extremely  low  permeability. 

•  The  NTS  experience  cannot  be  directly  translated  to 
proposed  deep  based  facilities  because  of  the  unique 
geologic  and  hydrologic  regime  at  NTS  and  because  of  the 
distinctive  characteristics  of  the  contained  burst 
loadings. 

—  The  Rainier  Mesa  tuff  with  its  perched  water  table, 
poor  hydraulic  interconnectivity  of  fractures, 
irregular  fracture  saturation  and  very  low  per¬ 
meability  is  not  characteristic  of  sites  being  con¬ 
sidered  for  deep  basing. 

—  There  are  important  differences  between  the  con¬ 
tained  burst  NTS  loadings  and  the  surface  or  shallow 
buried  burst  threat  to  deep  basing,  particularly  in 
the  mechanisms  resulting  in  long  term  residual  stress 
and  pore  pressure  buildups.  Despite  the  differences, 
however,  both  types  of  loadings  produce  residual 
stress  and  pore  pressure  buildups  which  will  produce 
pore  water  migration. 

•  Both  sources  of  explosion  generated  water  observed  at  NTS 
are  a  serious  threat  to  deep  based  systems  at  candidate 
sites  exemplified  by  the  Generic  Mountain  C  hydrologic/ 
geologic  setting. 

—  At  a  site  where  all  fractures  are  generally 
saturated,  explosion  mobilized  fracture  flow  could 
be  a  more  general  and  serious  problem  than  at  NTS. 
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However,  there  are  little  or  no  data  on  properties 
controlling  fracture  flow  at  Mountain  C,  such  as 
fracture  storage  capacity,  hydraulic  inter¬ 
connectivity,  and  the  mechanical  and  flow  properties 
of  the  fractures. 

—  Simplified  two  phase  calculations  of  flow  produced 
by  explosively  generated  residual  stresses  and  pore 
pressures  indicate  that  such  flow  is  a  very  serious 
potential  threat  to  deep  based  systems  in  the  Generic 
Mountain  C  setting.  Even  through  pore  water  mi¬ 
gration  was  not  a  serious  threat  at  NTS,  a  similar 
residual  stress  loading  of  the  much  more  permeable 
Generic  Mountain  C  geology  can  produce  flow  volumes 
orders  of  magnitude  greater  than  those  experienced 
at  NTS. 


RECOMMENDATIONS 

•  The  magnitude  of  the  threat  of  explosively  generated 
flooding  of  deep  based  systems  should  be  more 
realistically  defined  at  generic  sites  of  interest. 
Better  threat  definition  will  require  analysis  of 

both  the  fracture  water  and  pore  water  migration  threat 
mechanisms.  It  should  also  include  analysis  of  alter¬ 
natives  for  eliminating,  reducing  and/or  handling  such 
f looding. 

•  Reliable  predictions  of  fracture  flow  require  the 
development  of  analytic  tools  and  most  importantly  data 
describing  the  geologic/hydrologic  properties  and 
response  parameters  governing  such  flow. 
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—  Acquisition  of  this  descriptive  data  will  require 
both  laboratory  and  field  programs  designed  to 
measure  simulated  and  actual  dynamic  and  residual 
pore  pressure  response  and  flow,  and  the  develop¬ 
ment  of  programs  and  tests  to  measure  in  situ 
fracture  properties  governing  fracture  flow. 

—  Development  of  analytic  tools  should  include  two 

phase  codes  capable  of  modeling  the  dynamic  response 
of  saturated  fractures  and  the  resultant  residual 
pore  pressure  buildups  and  flows.  These  codes  should 
initially  be  used  to  conduct  parameter  studies  in 
which  fracture  flow  parameters  are  varied  over  the 
entire  range  of  interest  to  deep  based  systems. 

•  Additional  work. to  more  realistically  define  the  threat 
from  pore  water  migration  should  be  mostly  analytic, 
since  the  governing  material  parameters  are  much  better 
defined  than  those  controlling  fracture  flow.  It  should 
be  noted  that  nearly  all  the  refinements  suggested  below 
will  tend  to  reduce  the  flows  from  the  preliminary 
calculations  performed  in  this  initial  study.  Additional 
analytic  calculations  using  both  one  and  two  dimensional 
two  phase  codes  are  required.  One  dimensional  cal¬ 
culations  should  be  used  to  study: 

—  the  use  of  pressure  grouting  to  reduce  pore  water 
migration ; 

—  more  realistic  material  models,  including  the 
influence  of  shock  conditioning; 
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—  the  influence  of  partial  saturation  (as  reported 
for  Generic  Mountain  C)  on  flow; 

—  the  influence  of  limited  aquifer  boundaries  on  flow 
and  flow  volumes. 

Two  phase  calculations  using  currently  available  codes 

should  be  used  to: 

—  perform  more  realistic  modeling  of  the  ONETON 
experiment; 

—  realistically  model  the  response  of  depressed 
water  tables  (phreatic  surfaces)  in  the 
vicinity  of  tunnels; 

—  simulate  more  realistic  nonsymmetric  loading 
conditions  from  threat  surface  and  shallow  buried 
bursts  on  aquifers  having  locations,  properties  and 
boundary  conditions  representative  of  actual  field 
conditions. 
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SECTION  2 


GEOLOGIC  AND  HYDROLOGIC  SETTING,  RAINIER  MESA 


GEOLOGIC  SETTING 

Rainier  Mesa,  shown  in  Figure  2.1,  is  located  in  the 
north  central  region  of  the  Nevada  Test  Site.  The  mesa  is 
centered  at  about  37®  12'  N  latitude  and  116®  13'  W  longitude;  or 
at  about  890,000  ft  N  and  630,000  ft  E  in  the  center  zone  of  the 
Nevada  state  coordinate  system.  Rainier  Mesa  is  overlain  by 
volcanic  cap  rock  with  a  surface  elevation  between  7400  and  7680 
ft  above  sea  level.  The  top  of  the  mesa  is  more  than  2500  ft 
above  Yucca  Flat  in  the  nearby  basin.  The  mesa  runs  north-south. 
Width  of  the  cap  rock  is  approximately  1.5  miles,  length  about  3 
miles  and  area  about  4.4  square  miles.  According  to  Thordarson 
(1965),  the  average  rainfall  on  the  mesa  was  7.5  in/yr  in  the  5 
year  period  between  1959  and  1964. 

Rainier  Mesa  is  a  remnant  of  an  eroded  volcanic  plateau. 
The  volcanic  tuffs  making  up  the  mesa  are  of  moderately  recent 
(Miocene)  to  very  recent  (Pliocene)  origin.  Thordarson  (1965) 
identifies  11  layered  tuff  units  atop  the  underlying  sedimentary 
and/or  metamorphic  rocks.  The  tuff  layers,  2000  to  5000  ft  in 
thickness,  are  relatively  flat,  with  dips  of  generally  less  than 
25®.  These  overlie  much  older  sedimentary/metamorphic  dolomites, 
argillites  and  quartzites  of  Paleozoic  to  late  Precambrian  age. 
Prior  to  deposition  of  the  tuff,  the  underlying  rocks  had 
undergone  extensive  deformation  and  are  highly  fractured,  folded 
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and  faulted.  The  volcanic  tuff  was  deposited  atop  the  eroded 
surface  of  these  older  rocks;  hence,  the  tuff  beds  tend  to 
parallel  the  predepositional  topography.  The  tuffs  tended  to 
fill  in  the  underlying  valleys  so  that  the  bedding  planes  in  the 
later  deposits  are  flatter  than  those  in  the  earlier  ones. 

Three  major  types  of  tuff  make  up  Rainier  Mesa,  zeolitic- 
bedded  tuff,  fr iable-bedded  tuff  and  partially  welded  to  welded 
tuff.  The  Rainier  tuff  was  originally  composed  primarily  of  fine 
grained  pumice  and  glass  shards.  Tuff  is  formed  by  the  release 
of  gases  in  molten  lava  resulting  in  expansion  and  frothing  of 
the  lava.  As  the  froth  breaks  to  the  surface,  the  pumice  and 
glass  fragments  are  sometimes  ejected  high  into  the  air  by 
volcanic  explosions  and  deposited  in  the  form  of  ash  fall  tuff. 
Alternatively,  froth  may  foam  down  the  side  of  a  volcano  in  a 
glowing  avalanche  and  be  deposited  as  ash  flow  tuff.  Under 
certain  conditions  the  tuff  welds  itself  together  after 
deposition.  Depending  on  the  weight  of  overburden  and  degree  of 
melt  in  _the  freshly  deposited  material,  various  degrees  of  welded 
or  partially  welded  tuff  can  be  formed  during  cooling. 

The  tuff  units  identified  by  Thordarson  (1965)  making  up 
Rainier  Mesa  are  listed  in  Table  2.1.  The  underground  weapons 
effects  tests  are  conducted  in  the  lowermost  or  tunnel  bed  tuffs, 
members  Tl  through  T4.  These  overlie  the  older  sedimentary  rocks 
and  are  made  up  of  zeolitic  bedded  tuffs.  In  addition  to  the 
tunnel  bed  tuffs,  the  overlying  lower  portion  of  the  Grouse 
Canyon  member  of  the  Indian  Trail  formation  and  the  lower  portion 
of  the  Paintbrush  formation  are  also  zeolitic  tuffs.  Thordarson 
(1965)  lists  the  aggregate  thickness  of  these  zeolitic  units  at 
between  800  and  1200  ft.  The  pumice  and  glass  particles  of  the 
original  ash  flow  tuffs  in  these  units  have  been  highly  altered 
to  form  various  zeolitic  minerals.  These  minerals  form  a  nearly 
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impermeable  matrix  around  the  remaining  non-zeolitic  constituents 

such  as  quartz,  feldspar,  biotite,  etc.  The  resultant  zeolitic 

tuff  is  a  highly  porous  but  very  impermeable  material.  As  shown 

in  Table  2,1,  Thordarson  (1965)  gives  average  interstitial 

porosities  of  the  zeolitic  tunnel  bed  tuffs  from  the  Ul2e  tunnel 

complex  of  from  25  to  38%.  On  the  other  hand,  average 

interstitial  permeabilities  range  form  9.4  x  10“^  to  1.9  x 
-•8 

10  cm/s,  in  the  very  low  to  practically  impermeable  range. 

There  appears  to  be  no  relationship  between  porosity, 
permeability  and  grain  size  in  the  zeolitic  tuff. 

The  second  major  tuff  category  forming  Rainier  Mesa  is 
the  friable  bedded  tuff.  As  shown  in  Table  2,1,  this  tuff  occurs 
in  the  lower  portion  of  the  Grouse  Canyon  member  and  comprises 
the  bulk  of  the  Paintbrush  Tuff  formation.  The  friable  tuff  is  a 
porous,  weak,  ash  fall  tuff  in  which  the  pumice  and  glass 
shards  are  largely  unaltered.  The  porosity  tends  to  average 
somewhat  higher  than  the  zeolitic  tuff;  samples  from  the  U12b 
tunnel  complex  averaged  40%  porosity.  The  average  interstitial 
permeability  of  the  friable  tuff  is  much  larger  than  that  in  the 
zeolitic  tuff.  Permeability  measurements  from  Emerick  and  Houser 
(1962)  average  1.75  x  10”^  cm/s;  lower  values  are  given  by 
Bowers  (1963)  with  an  average  value  of  2.4  x  10”^  cm/s.  The 
latter  measurements  were  made  using  air  rather  than  water.  Data 
from  Keller  (1960)  indicate  that  permeability  to  air  is  2  to  2Ci 
times  higher  than  the  permeability  to  water.  The  above 
permeabilities  fall  in  the  low  permeability  range  (10”^  to 
10  ^  cm/s)  according  to  Lambe  and  Whitman  (1969). 

The  welded  and  partially  welded  tuffs  are  among  the  least 
porous  and  least  permeable  rocks  within  Rainier  Mesa.  As  shown 
in  Table  2.1,  these  tuffs  comprise  the  Rainier  Mesa  cap  rock,  the 
Tiva  Canyon  and  part  of  the  Stockade  Wash  members  of  the 
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Paintbrush  formation  and  the  upper  Grouse  Canyon  and  Tub  Spring 
members  of  the  Indian  Trail  formation.  The  welded  tuffs  were 
deposited  as  incandescent  ash  falls  which  were  subsequently 
self-welded  under  their  own  weight  and  heat.  Cooling  fractures 
are  abundant  in  the  welded  tuff,  with  the  highest  density  of 
fractures  in  the  denser  more  highly  welded  tuffs.  According  to 
Table  2.1,  porosity  of  the  Rainier  cap  rock  averaged  only  14% 
while  that  in  the  Grouse  Canyon  welded  tuff  averaged  19%. 
Interstitial  permeability  in  the  cap  rock  averaged  4.7  x  10“^ 
cm/s  and  the  one  sample  of  Grouse  Canyon  welded  tuff  reported  had 
a  permeability  of  2.8  x  10  °  cm/s.  These  permeabilities  are 
in  the  very  low  to  practically  impermeable  range  described  by 
Lambe  and  Whitman  (1969). 

Beginning  in  moderately  recent  Miocene  time  and 
continuing  into  the  very  recent  Pliocene,  the  Rainier  Mesa  region 
has  been  subjected  to  normal  faulting  resulting  in  the  depression 
of  the  adjacent  Yucca  Valley  to  the  east  and  Fortymile  Canyon  to 
the  west.  In  conjunction  with  this  faulting  the  Rainier  Mesa 
tuff  has  undergone  extensive  normal  faulting  and  joint 
development.  Faults  are  very  steep  to  nearly  vertical  and 
exhibit  varying  amounts  of  normal  displacement.  Most  faults  have 
stratigraphic  normal  displacements  of  inches,  though  major  faults 
show  displacements  to  tens  of  feet.  Major  faults  transect  all 
the  tuff  beds  within  the  mesa,  including  the  cap  rock  at  the 
surface.  Minor  faults  are  much  more  numerous  than  major  faults, 
with  many  extending  less  than  300  ft.  Joint  and  fault  spacing 
vary  dramatically  within  relatively  small  local  areas,  from  tens 
of  feet  to  inches.  Openings  in  faults  and  joints  also  vary 
considerably.  Fault  openings  of  up  to  6  in  are  reported  by 
Thordarson  (1965),  but  in  most  instances  faults  are  relatively 
tight  and  may  be  filled  with  clay  gouge.  Joints  are  also 
generally  closed,  but  openings  of  up  to  several  inches  are 
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observed.  Thordarson  also  observes  that  fracture  openings  vary 
irregularly  throughout  the  rock  masSf  with  fractures  that  are  open 
several  inches  at  one  point  being  tightly  closedonly  a  few  feet 
away . 


HYDROLOGIC  REGIME 

Ground  water  within  Rainier  Mesa  is  concentrated  in  the 
zeolitic  tunnel  beds  1  through  4  of  Table  2.1.  The  very  low 
permeability  of  the  tunnel  bed  tuffs  traps  water  which  originated 
as  rain  on  the  mesa  surface  and  which  percolates  downward, 
primarily  through  fractures  in  the  overlying  tuffs.  The  water 
trapped  in  the  zeolitic  tuffs  forms  a  perched  water  table 
approximately  2000  ft  above  the  true  regional  water  table  which 
lies  well  beneath  the  floor  of  the  surrounding  valleys. 

Thordarson  (1965)  classes  the  zeolitic  bedded  tuff  as  a  fractured 
aquitard.  There  is  no  appreciable  flow  of  water  through  the 
interstices  of  the  zeolitic  tuff.  Migration  of  water  downward 
through  the  tunnel  bed  layers  is  predominantly  through  the 
generally  tight  and  poorly  connected  fractures. 

The  top  of  the  perched  water  table  shows  considerable 
variation  in  elevation  due  to  the  poor  hydraulic  inter¬ 
connectivity  of  the  fractures.  It  lies  in  the  top  of  the 
zeolitic  tunnel  beds,  generally  within  a  few  hundred  feet  of  the 
6000  ft  elevation.  Figure  2.2  shows  the  location  of  various 
wells  from  which  Thordarson  (1965)'  inferred  the  elevation  of  the 
zone  of  fracture  saturation  in  the  vicinity  of  the  E  tunnel 
complex  (Ul2e).  The  wells,  shafts,  and  exploration  holes 
discussed  by  Thordarson  are  listed  in  Table  2.2.  Only  two  of  the 
holes,  Hagestad  1  and  Ul2e.06-1,  were  drilled  from  the  top  of 


13 


the  mesa.  The  others  were  drilled  or  excavated  from  within  the  E 
or  B  tunnel  complexes.  In  addition,  in  no  instance  was  the  depth 
to  the  top  of  the  fracture  saturation  surface  measured  in  an 
unambiguous  manner.  Of  the  two  holes  from  the  mesa  surface, 
U12e.06-1  was  evidently  drilled  through  the  water  bearing 
zeolitized  tuff  beds  into  the  underlying  unsaturated  dolomite 
without  measuring  the  perched  water  table  in  the  tunnel  bed 
tuffs.  The  Hagestad  1  test  hole  was  cased  to  its  full  depth  of 
1941  ft  and  cemented.  The  casing  was  then  perforated  at  selected 
depths  allowing  water  to  fill  the  hole.  The  depth  of  water  in 
the  hole  was  a  function  of  the  location  of  the  casing 
perforations.  The  highest  elevation  of  water  in  the  hole  was 
6039  ft,  but  Thordarson  admits  that  perforations  at  other 
elevations  might  have  filled  the  hole  further. 

The  remaining  holes  and  shafts  listed  in  Table  2.2  were 
excavated  or  drilled  from  within  the  E  and  B  tunnel  complexes 
into  the  underlying  tuffs.  In  the  case  of  holes  U12e.03-1  and 
Ul2e.M~l,  water  poured  from  the  hole  up  into  the  tunnel.  These 
holes  were  plugged  off  and  the  elevation  of  the  resulting  head 
computed  from  the  pressure  in  the  hole.  In  all  measurements  made 
from  within  the  tunnels,  the  very  presence  of  the  tunnel  had 
influenced  the  water  table  elevation  in  that  area  as  evidenced  by 
the  large  volumes  of  water  which  were  produced  during  excavation 
(production  of  water  is  discussed  in  Section  3).  For  the  above 
reasons,  the  water  table  elevations  of  Table  2.2  should  be 
considered  minimum  elevations.  Elevations  prior  to  tunnel 
construction  and  elevations  in  the  vicinity  of  the  cased  drill 
holes  were  probably  somewhat  higher  than  indicated  in  the  table. 
At  any  rate,  in  the  vicinity  of  E  Tunnel  the  elevation  of  the 
zone  of  fracture  saturation  appears  to  be  somewhat  in  excess  of 
6000  ft,  with  considerable  variation  within  the  range  of  at  least 
6000  ft  to  6200  ft  elevation. 
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Thordarson  (1965)  attributes  the  large  variability  in 
water  table  to  the  variations  in  fracture  permeability  and 
hydraulic  interconnectivity  within  the  zeolitized  tuff.  In 
areas  where  fracture  frequency  and  permeability  are  relatively 
high,  one  might  expect  downward  seepage  of  ground  water  to  be 
enhanced  and  the  water  table  to  be  somewhat  depressed  relative  to 
more  impermeable  adjacent  areas.  Winograd  and  Thordarson  (1975) 
present  the  simplified  cross  section  shown  in  Figure  2.3  to 
schematically  illustrate  the  perched  water  table  in  the  vicinity 
of  E  tunnel.  Beneath  the  eastern  two  thirds  of  the  mesa,  the 
underlying  dolomite  (lower  carbonate  aquifer)  is  unsaturated  and 
lies  well  above  the  regional  water  table.  The  water  within  the 
tunnel  bed  tuffs  (tuff  aquitard)  is  perched  above  the  unsaturated 
dolomite.  The  variable  water  table  in  the  region  of  E  Tunnel  is 
represented  in  the  triangular  shaped  fracture  zones  1  and  4.  At 
location  1  variable  levels  of  perched  water  in  the  fracture  zones 
are  represented  by  the  variable  black  shading  of  the  schematic 
fracture  zones.  In  this  region,  the  water  table  lies  even  with, 
or  below,  the  level  of  E  Tunnel.  At  location  4,  the  perched  water 
table  lies  well  above  the  tunnel  level. 

Beneath  the  western  third  of  the  mesa  the  underlying 
dolomite  dips  beneath  the  regional  water  table  and  the  tuff 
aquitard  dips  to  the  level  of  the  regional  water  table.  This 
creates  a  continuous  saturated  zone  above  the  regional  water 
table  which  is  termed  a  semiperched  zone.  The  semiperched  water 
table  extends  to  the  base  of  the  overlying  welded  tuff  aquifer 
which  dips  below  the  6000  ft  elevation  in  this  region. 

Immediately  south  of  the  mesa  at  the  location  of  Well 
87-62  (Test  Well  1  in  Figure  2.2),  the  top  of  the  tuff  aquitard  is 
at  an  elevation  of  5931  ft,  225  ft  beneath  the  well  head 
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elevation  of  6156  ft.  The  top  of  the  perched  water  table  is  at 
an  elevation  of  5746  ft/  only  410  ft  below  the  well  head  and  185 
ft  beneath  the  top  of  the  aquitard.  The  perched  water  table  at 
this  location  is  1560  ft  above  the  regional  water  table#  which  has 
an  elevation  of  4186  ft  (Winograd  and  Thordarson#  1975). 

Thordarson  (1965)  notes  that  these  water  levels  from  Test  Well  1 
are  the  only  true  unambiguous  measurements  of  the  static  water 
level  in  the  region.  It  is  interesting  to  note  that  if  the  top 
of  the  perched  water  table  were  also  about  200  ft  beneath  the  top 
of  the  tuff  aquitard  in  the  vicinity  of  E  Tunnel,  it  would  be 
approximately  150  ft  above  the  tunnel  elevation. 

Thordarson  (1965)  maintains  that  the  tunnel  bed  tuffs  are 
fully  saturated  intersti tially ,  not  only  within  the  zone  of 
fracture  saturation,  but  also  hundreds  of  feet  above  it  as  well. 

He  notes  seeps  of  water  which  were  found  in  the  Ul2b  tunnel 
system  at  an  elevation  of  6600  ft,  some  400  to  600  ft  above  the 
level  of  fracture  saturation  in  this  region.  Extensive  analysis 
of  tunnel  bed  tuffs  since  that  time  have  indicated  that  there  is 
some  variability  in  the  degree  of  interstitial  saturation.  In 
some  locations,  the  tunnel  bed  tuffs  appear  to  contain  a  few 
percent  air  voids,  though  specification  of  the  exact  location  of, 
and  amount  of,  these  air  voids  has  been  very  illusive. 

Based  on  an  estimated  annual  ground  water  recharge  for 
Rainier  Mesa  from  the  measured  annual  rainfall,  and  on  measured 
hydraulic  gradients  from  Test  Well  1  and  Hagestad  1,  Thordarson 
(1965)  estimates  an  average  vertical  permeability  for  the  tuff 
aquitard  in  the  9.43  x  10”®  to  2.36  x  10"®  cm/sec  range. 

Note  that  this  range  agrees  well  with  the  interstitial 
permeabilities  for  the  tunnel  bed  tuffs  reported  in  Table  2.1. 
Since  it  agrees  with  interstitial  permeabilities  from  laboratory 
samples  it  would  appear  that  Thordarson 's  hypothesis  that  the 
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effective  permeability  is  governed  by  the  fractures  is  incorrect 
However,  experience  in  mining  the  tunnels  beneath  Rainier  Mesa, 
summarized  in  the  following  section,  shows  that  the  effective 
permeability  is  indeed  governed  by  the  fracture  flow.  Several 
factors  may  account  for  the  apparent  inconsistency  between  the 
laboratory  permeabilities  and  Thordarson's  derived  in  situ 
permeabilities.  These  factors  are: 

•  The  average  permeabilities  presented  in  Table 

2.1  may  not  be  representative  of  the  tunnel  bed  tuffs; 

•Thordarson's  assumed  portion  of  the  annual  rainfall 
going  into  ground  water  recharge  is  too  low  (this  was 
based  on  data  from  other  sites); 

•  the  measured  hydraulic  gradients  upon  which  Thordarson's 
analysis  was  based  are  not  representative  of  those  in 
the  tuff; 

•  or  some  combination  of  these  factors. 

Measurement  of  in  situ  or  effective  permeabilities  in  the  tunnel 
bed  tuffs  is  a  subject  needing  further  investigation. 
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Table  2.1.  Rainier  Mesa  lithology  and  properties  (after  Thordar  965) 


Table  2.2.  Data  on  elevation  of  fracture  satxiration  in 
the  vicinity  of  E  Tunnel  (from  Thordarson, 
1965)  . 


Hole 

Total  Depth 

Hole  Elevation 

Inferred  Elev, 
of  Fracture 
Saturation 

U12e03-1 

834  ft  deep 

El  6150  ft 

6167  ft 

U12e.M-l 

1501  ft  deep 

El  6158  ft 

6184  ft 

lagestad  1 

1941  ft  deep 

El  7485  ft 

6039  ft 

U12e.06-1 

3114  ft  deep 

El  7573  ft 

<  4643  ft 

Shaft  in  U12e.07 

not  reported 

not  reported 

6033  ft 

Shaft  in  U12b.07 

not  reported 

not  reported 

6147  ft 

19 


1^3.000  rf( 


Figure  2 


hii'  ' 

W«l4l«d*luff  •HWifcf 


E3 

Bc^tfcd'lyff  ytiuifcf 


Tuff  m^u%iu4 


Lower  ctflHNicie  Mfuifer 


f-XTLANATION 

T 

Tnt  Hole 

FRACTUKK  TYfKS 

^  hm  ocMIoii  lOiAlt  or  Uulu.  sooir  cofiuin  wuitr.  ofhert 
eoi^y 

^  Foviu.  clofod  •hove  luff  •i|Miuril*€»rhonaif  aiiuifvr  con- 
laci.  fully  and  partly  taluralad 
Q)  Fault,  opafi  at  hoitom,  amply 

^  Fault,  wiih  pmch-and-awall  atruciura.  partly  ■aCuraiad 
Parahad  and  lamiparchcd  water  »  Nack 


3.  Schematic  section  view  of  perched  and  semiperched 
groundwater  in  the  vicinity  of  E  Tunnel,  Rainier 
Mesa  (from  Winograd  and  Thordarson,  1975). 
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SECTION  3 


OCCURRENCE  OF  WATER  DURING  AND  AFTER  MINING 


OVERVIEW 

Flow  of  water  into  the  tunnel  complexes  in  the  zeolitized 
tuffs  beneath  Rainier  Mesa  has  been  documented  in  a  piecemeal 
fashion  for  many  years.  A  brief  look  at  some  of  the  descriptive 
and  quantitative  features  of  this  flow  is  very  helpful  in 
understanding  the  occurrence  and  migration  of  the  water  in  the 
tunnel  bed  tuffs.  It  is  also  very  helpful  in  formulating 
hypotheses  whicn  explain  the  tunnel  flooding  which  sometimes  has 
occurred  adjacent  to  underground  nuclear  and  HE  detonations. 

A  plan  view  of  the  tunnel  complexes  beneath  Rainier  Mesa 
is  shown  in  Figure  3.1.  The  major  complexes  are  the  F,  E,  B,  N 
and  T  tunnel  complexes;  except  for  B  Tunnel  all  were  mined  in  the 
zeolitized  tunnel  bed  tuff  units  1  through  4  of  Table  2.1.  B 
Tunnel,  with  a  portal  elevation  of  6606  ft,  was  mined  in  the 
overlying  Lower  Grouse  Canyon  member  of  the  Paintbrush  Tuff. 
Significant  flows  of  water  from  faults  and  joints  were 
encountered  in  the  E,  N,  and  T  complexes  during  mining.  In  some 
instances  flow  continued  well  after  mining  operations  had  moved 
elsewhere.  In  G  Tunnel,  toward  the  extreme  south  end  of  Rainier 
Mesa,  only  a  few  small  seeps  were  encountered  during  the  mining. 
In  E  Tunnel,  further  into  the  mesa,  but  having  the  same  portal 
elevation  (6115  ft)  as  G  Tunnel,  numerous  seeps  and  many  moderate 
flows  were  encountered  during  mining.  Only  a  few  small  seeps 
were  encountered  in  B  Tunnel  which  was  mined  in  the  relatively 
permeable  Paintbrush  Tuff  above  the  perched  water  table.  In  N 
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and  T  Tunnels,  with  portal  elevations  of  6024  ft  and  5600  ft 
respectively,  numerous  small  seeps  and  many  moderate  flows  were 
encountered. 

The  remainder  of  this  section  will  highlight  features  of 
the  flow  encountered  in  the  E,  N  and  T  complexes  as  well  as 
observations  related  to  the  hydraulic  interconnectivity  and 
capacity  of  the  fractures. 


E  TUNNEL 

Thordarson  (1965)  and  Clebsch  (1960)  provide  excellent 
descriptions  of  the  seeps  and  flows  encountered  during  the  mining 
of  E  Tunnel.  From  August  1958  through  December  1963  an  estimated 
30  to  50  million  gallons  (90  -  150  acre  ft)  of  water  flowed  from 
the  E  Tunnel  portal.  This  water  was  discharged  from  the  faults 
and  joints  exposed  during  the  mining  operations.  A  detailed 
layout  of  the  E  Tunnel  complex  is  shown  in  Figure  3.2.  IncjLuded 
on  this  map  are  the  locations  of  faults  intersecting  the  various 
drifts.  A  total  of  about  110  faults  and  5000  joints  were  mapped 
during  the  mining  of  approximately  19,000  ft  of  drifts  in  E 
Tunnel.  Of  these,  about  50%  to  60%  of  the  faults  produced  most 
of  the  water,  while  about  2%  of  the  joints  yielded  a  minor 
portion  of  the  water. 

Of  the  total  30  to  50  million  gallon  flow  from  E  Tunnel, 
the  largest  contribution  probably  came  from  the  U12e.02  LOGAN 
drift  and  the  Ul2e.05  BLANCA  drift.  Between  the  middle  of  August 
and  end  of  October  1958,  approximately  10  million  gallons  of 
water  flowed  from  these  two  drifts  alone.  During  the  ensuing  3 
year  period  ending  in  December  1961,  about  15  to  30  million 
additional  gallons  of  water  flowed  out  of  the  E  Tunnel  complex  as 


24 


the  main  drift  was  lengthened  and  adjacent  drifts  were  added. 
During  the  final  period  reported  (Thordarson,  1965),  between 
December  1961  and  December  1963,  another  5  to  10  million  gallons 
was  measured,  most  of  this  from  mining  of  the  U12e.03  and  U12e.06 
drifts. 


Thordarson  (1965)  reports  that  about  half  the  occurrences 
of  water  were  from  in  or  near  faults,  and  that  most  instances  of 
flows  of  more  than  5  gallons  per  minute  (gpm)  were  directly  out 
of  faults.  The  larger  flows  of  fracture  water  from  the  faults  in 
E  Tunnel  ranged  from  5  gpm  to  20  gpm.  Typically,  after  a 
fracture  was  penetrated  by  a  tunnel  there  was  a  maximum  initial 
discharge  which  decreased  gradually  to  a  small  seep  within  a  few 
days.  Within  a  few  weeks  or  months  most  of  the  fractures  had 
drained  completely;  however,  Clebsch  {I960)  reports  that  water 
dripped  from  some  fractures  for  a  period  of  2  years  or  more. 

The  rapid  decrease  in  flow  is  indicative  of  poorly 
connected  fracture  surfaces  within  the  tunnel  bed  tuffs.  Both 
the  fact  that  the  fractures  drained  relatively  quickly  and  that 
closely  spaced  fractures  flowed  strongly  following  drainage  of 
adjacent  fractures  attest  to  this  conclusion.  There  also 
appeared  to  be  a  correlation  between  flow  and  the  number  or 
density  of  fractures.  Thordarson  (1965)  notes  that  the  driest 
drifts  in  the  E  Tunnel  complex  were  the  Ul2e.01  and  (Jl2e.07 
drifts  (see  Figure  3.2)  which  also  contained  the  fewest 
fractures. 

The  variability  in  flow  and  .correlation  with  fracture 
density  is  illustrated  in  the  plot  of  flow  vs.  time  for  the 
Ul2e.05  BLANCA  and  Ul2e.02  LOGAN  drifts  shown  in  Figure  3.3.  The 
data  record  the  total  flow  out  of  both  drifts  as  reported  by 
Clebsch  (1960)  during  the  mining  which  progressed  as  shown  in  the 
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bottom  portion  of  the  figure.  Flow  varied  from  a  low  of  about  30 
gpm  to  a  high  of  about  500  gpra  during  the  mining,  with 
considerable  variation  as  flow  volumes  declined  in  exposed 
fractures  and  were  suddenly  increased  as  new  fractures  were 
exposed.  Within  two  weeks  after  the  completion  of  the  drifts  the 
total  flow  had  decreased  to  about  20  gpm.  The  greatest  flow 
occurred  during  penetration  of  intensely  jointed  and  faulted  tuff 
of  Tunnel  Bed  3  in  the  Ul2e.05  BLANCA  drift.  Subsequent  high 
discharge  values  were  measured  during  penetration  of  fractures  in 
the  Ul2e.02  LOGAN  drift.  At  the  time  of  the  last  measurement 
most  of  the  flow  was  from  the  LOGAN  drift. 

Thordarson  (1965)  notes  that  there  were  tens  of  open 
fractures  within  the  E  Tunnel  complex  which  contained  no  water, 
despite  the  fact  that  fractures  on  either  side  were  full  of 
water.  The  presence  of  these  empty  fractures  is  explained  as 
either  due  to  their  complete  isolation  from  the  surrounding 
saturated  fractures  which  receive  recharge,  or  due  to  their  being 
open  below  and  connected  through  discharge  channels  to  the 
underlying  regional  water  table. 

He  also  notes  that  locally  there  are  certain  joint 
orientations  which  contain  all  of  the  fracture  water.  For 
instance,  at  a  range  of  3000  to  3500  from  the  portal  in  the  main 
drift,  water  occurred  only  in  the  NE-SW  striking  principal  joint 
set.  No  water  flowed  from  the  NW-SE  joints  forming  the  minor  set. 
In  contrast,  water  occurred  in  the  NW-SE  major  joint  set  in  the 
Ul2e.03b  drift,  but  there  was  no  discharge  from  the  NE-SW  minor 
joint  set.  Such  contrasts  indicate  that  there  is  poor  hydraulic 
interconnectivity  between  the  major  and  minor  joint  sets  in  these 
areas. 
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Thordarson  concludes  that  production  of  water  in  the  U12e 
tunnel  complex  was  governed  by  the  extent,  density,  opening  width 
and  interconnectivity  of  the  steeply  dipping  faults  and  joints  in 
the  zeolitized  tuff  surrounding  the  tunnel.  Most  of  the'  faults  in 
the  area  appear  to  be  local,  perhaps  en  echelon  faults,  limited  in 
extent  to  less  than  100  to  300  ft.  For  example,  note  that  most  of 
the  densely  spaced  faults  along  the  Ul2e.05  BLANCA  drift  in 
Figure  3,2  don't  intersect  either  the  adjacent  Ul2e.03  or  U12e.05 
drifts.  Openness  of  the  faults  varies  from  up  to  6  in  to  closed 
and  nearly  sealed  by  fault  gouge.  Joints  are  generally  closed, 
but  in  some  areas  are  open  several  inches  at  one  point  and  are 
tightly  closed  within  just  a  few  feet. 

In  order  for  the  fractures  to  perch  the  downward 
migrating  fracture  water,  Thordarson  points  out  that  they  must 
be  closed  or  nearly  closed  at  some  underlying  location (s)  along 
their  strike.  Evidence  from  the  tunnels  suggests  that  some 
fractures  are  irregularly  open  and  pinched  shut  due  to  faulting 
action.  Other  fractures  appear  to  be  open  in  the  massive 
zeolitic  tuff,  but  closed  in  the  thin  underlying  beds  of  clayey 
tuff.  He  also  notes  that  clayey  gouge  or  other  fine  minerals 
could  be  deposited  by  the  downward  migrating  pore  water  to  seal 
off  open  faults. 

There  appears  to  be  no  evidence  that  interstitial  water 
contributes  to  flow  into  the  tunnels.  The  only  evidence  of  any 
interstitial  flow  are  localized  moist  spots  along  the  tunnel 
walls.  There  were  no  drips  or  other  flow  in  evidence.  These  wet 
spots  may  have  formed  due  to  seepage  controlled  by  particularly 
impermeable  zones  within  the  tuff.  In  these  locations,  the  rate 
of  evaporation  in  the  ventilated  tunnels  evidently  equalled  the 
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rate  of  water  production.  In  most  areas  the  tunnel  walls  were 
dry,  further  testimony  to  the  impermeability  of  the  zeolitic 
tuf  f . 


A  final  point  by  Thordarson  (1965)  concerning  flow  of 
interstitial  water  makes  an  interesting  contrast  to  the  data 
presented  from  N  Tunnel  in  the  following  subsection.  As  shown  in 
Figure  3.2,  there  is  a  syncline,  or  dip,  in  the  tunnel  beds  with 
its  trough  or  axis  running  NE-SW  across  the  E  tunnel  complex. 
Thordarson  notes  that  there  was  no  concentration  or  flow  of 
water  into  the  base  of  this  prominent  syncline  at  its 
intersection  with  any  of  the  three  drifts  it  crosses.  Only  minor 
seeps  were  observed  at  the  intersections  with  the  synclinal  axis 
and  these  were  from  fractures,  rather  than  due  to  flow  along 
bedding  planes  in  the  trough. 

N  TUNNEL 


Overall,  the  N  tunnel  complex  was  somewhat  drier  than  E 
Tunnel  during  and  after  mining.  However,  significant  water  and 
mining  problems  developed  during  the  mining  of  the  Ul2n.03  drift 
shown  in  Figure  3.4.  The  U12n.03  drift  was  designed  to  house  a 
line  of  site  pipe  for  an  underground  nuclear  test.  Mining  of  the 
drift  began  on  18  April  1966  and  was  completed  on  22  May  1967. 

The  drift  runs  N  26®W  from  the  Ul2n  extension,  a  distance  of  2166 
ft.  The  elevation  at  the  portal  end  of  the  drift  was  6067  ft  and 
the  drift  was  mined  upward  at  a  0.5%  grade.  The  entire  drift  is 
within  Tunnel  Bed  4  of  the  zeolitized  tuffs  shown  in  the  geologic 
section  of  Table  2.1.  Ege  et  al,  (1980)  describe  the  rock  around 
the  Ul2n.03  drift  as  an  ash  fall  tuff,  interbedded  with  reworked 
ash  fall  tuff  and  tuffaceous  sandstone.  The  tuff  is  bedded, 
zeolitized,  and  in  places  altered  to  a  high  clay  content.  The 
drift  crosses  the  Aqueduct  Syncline  approximately  1250  ft  from 
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the  main  drift.  The  Aqueduct  Syncline  is  a  major  syncline  with 
its  trough  oriented  nearly  perpendicular  to  the  Ul2n.03  drift. 


A  section  and  plan  view  of  the  U12n.03  drift  in  the 
vicinity  of  the  Aqueduct  Syncline  is  shown  in  Figure  3.5.  The 
axis  of  the  syncline  is  near  survey  station  13-*-00.  Bedded  tuff 
subunits  Tt4kr  Tt4J  and  Tt4H,  subunits  of  Tunnel  Bed  4,  are 
pictured  sloping  gently  upward  on  either  side  of  the  synclinal 
axis.  The  drift  is  within  the  Tt4k  subunit  for  approximately  320 
ft  SE  of  the  synclinal  axis  and  for  about  450  ft  on  the  NW  side  of 
the  axis. 


Ege  et  al.  (1980)  report  that  the  Tt4k  tuff  between 
stations  10+00  and  16+75  ft  has  been  strongly  altered  to  clay  by 
the  action  of  ground  water  which  has  collected  in  the  region. 

X-ray  analysis  indicated  the  presence  of  calcium-montmor illonite 
and  mechanical  analyses  showed  the  rock  to  be  highly  plastic  with  a 
low  unconfined  compressive  strength.  Unconfined  strengths  of 
core  samples  from  the  U12n.03  UG-3  drill  hole,  shown  in  Figure 
3.5,  averaged  1410  psi,  with  a  range  in  strengths  of  from  290  psi 
to  2300  psi  and  a  standard  deviation  of  670  psi.  Ege  et  al. 
estimate  that  maximum  stresses  in  the  walls  of  the  Ul2n.03  drift 
due  to  overburden  loads  and  stress  concentrations  are 
approximately  2200  psi.  As  a  result  of  the  concentrated 
overburden  stresses  exceeding  the  rock  strength,  severe 
construction  and  support  problems  developed  in  the  highly  plastic 
rock.  Between  survey  stations  12+00  and  15+25  swelling  and 
squeezing  ground  deformed  steel  sets  and  popped  the  plates  from 
rockbolts  used  to  stabilize  the  drift.  Within  a  two  week  period 
in  the  summer  of  1966  lateral  deformations  of  the  tunnel  walls  of 
up  to  36  in  were  monitored. 
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Also  contributing  to  the  instability  of  the  U12n.03  drift 
in  the  vicinity  of  the  synclinal  axis  were  a  number  of  normal 
faults  which  intersected  the  drift  at  a  shallow  angle  as  shown  in 
Figure  3.5.  The  high  density  of  faults  in  this  region,  and  the 
fact  that  they  intersected  the  drift  at  very  shallow  angles 
substantially  exacerbated  the  movement  of  the  squeezing  and 
swelling  ground. 

A  final  problem  encountered  on  the  NW  side  of  the 
synclinal  axis  was  heavy  water  flow.  Flow  locations  are  shown  in 
the  plan  view  of  Figure  3.5.  Minor  to  moderate  flows  of  water  of 
up  to  5  gpm  were  encountered  issuing  from  joints  and  faults  during 
mining  near  the  synclinal  axis  in  July  and  August  of  1966.  On  21 
September  1966  heavy  ground  water  flows  of  65  gpm  were  encountered 
from  the  face  of  the  drift  near  the  floor.  Ege  et  al.  (1980) 
report  that  after  ten  days  of  dewatering  this  flow  had  decreased 
to  40  gpm  and  that  after  50  days  the  flow  had  dropped  to  25  gpm. 
The  water  was  issuing  from  three  fractures  within  a  fracture  zone 
beneath  the  drift.  These  fractures  were  open  from  1  to  3  in  and 
dipped  steeply  toward  the  drift  portal.  A  sump  was  built  in  the 
fracture  zone  and  water  was  pumped  from  the  zone  at  rates  varying 
between  40  and  60  gpm.  In  January  of  1968,  more  than  two  years 
after  excavation,  the  flow  rate  was  measured  at  8  gpm. 

Due  to  the  severe  water  problems  and  unstable  ground,  the 
Ul2n.03  drift  was  abandoned.  A  bulkhead  was  constructed  near  the 
drift  portal  and  for  many  years  the  drift  was  used  for  a  water 
supply.  In  1979  the  drift  was  reentered  for  a  limited  distance 
at  both  ends  in  conjunction  with  site  investigations  for  the 
MINERS  IRON  event.  In  April  of  that  year  total  flow  rate  from 
the  Ul2n.03  drift  was  0.35  gpm. 
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T  TUNNEL 


An  interesting  study  of  water  flow  and  pressure  buildup 
was  made  in  an  array  of  six  horizontal  drill  holes  running  north 
and  west  off  the  end  of  the  Ul2t  main  drift.  This  array 
consisted  of  NX  exploration  holes  ranging  in  length  from  1500  to 
3700  ft.  The  hole  layouts  are  shown  in  the  expanded  view  of  the 
NW  end  of  the  T  tunnel  complex  in  Figure  3.6.  All  holes  except 
U12t.03  UG-1  were  drilled  from  near  the  end  of  the  U12t  main 
drift.  Hole  Ul2t.03  UG-1  was  drilled  from  the  end  of  the  U12t.01 
Bypass  drift  on  a  bearing  of  approximately  N  65 ®W.  This  was  the 
longest  of  the  6  holes  and  the  first  to  be  drilled.  It  was 
drilled  throughout  most  of  1972,  being  completed  to  a  length  of 
3690  ft  in  early  November.  The  second  hole  to  be  completed  was 
Ul2t.03  UG-2  completed  to  a  length  of  1504  ft  in  December  1972. 
The  three  holes  sharing  a  common  terminus  near  the  end  of  the 
main  drift,  U12t.04  UG-1,  U12t.05  UG-1  and  Ul2t.06  UG-1,  were 
completed  in  that  order  in  May,  August  and  October  of  1973. 
Finally  hole  Ul2t.03  UG-3  was  completed  in  July  of  1974  on  a 
bearing  of  N  25®W  from  the  the  terminus  of  U12t  main  drift. 

The  HUSKY  PUP  line  of  site  drift  (U12t.03)  was  rained  along  the 
path  of  the  U12t.03  UG-3  exploration  hole.  Predictions  of  ground 
water  flow  into  the  LOS  drift  during  mining  were  made  by  Hoover 
(1974)  based  on  the  ground  water  inflows  into  the  U12t.03  UG-3 
exploration  hole.  An  overall  description  of  water  production  in 
all  the  Ul2t  tunnel  drill  holes  is  given  in  a  memorandum  by 
Hoover  (1975). 

A  summary  of  the  peak  flows  and  pressures  measured  in  the 
6  U12t  drill  holes  is  shown  in  Table  3.1.  In  all  but  two  of  the 
holes,  peak  flow  occurred  after  completion  of  the  drilling. 
Examination  of  the  flow  logs  in  Hoover  (1975)  indicates  that 
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flow  originated  in  fracture  zones  within  each  hole.  The  drilling 
and  flow  logs  from  hole  U12t.04  OG-1  are  plotted  in  Figure  3.7  to 
illustrate  the  variability  in  flow  with  hole  length.  Water  was 
first  encountered  at  a  hole  length  of  320  ft  on  31  April  1973.  The 
cumulative  flow  built  slowly  to  10  gpm  on  17  May  at  a  length  of 
1209  ft.  Flow  suddenly  increased  to  25  gpm  as  the  hole  length 
was  extended  from  1209  ft  to  1245  ft  on  the  18th,  the  increase 
being  attributed  to  water  bearing  fractures  between  these  ranges. 
Flow  then  remained  steady  until  the  21st,  when  a  major  water  bearing 
fracture  zone  was  penetrated  between  1280  and  1400  ft.  Over  this 
interval  the  total  flow  increased  to  200  gpm.  Total  flow  then 
decreased  over  the  next  week  as  the  water  in  this  zone  was 
depleted.  Near  the  end  of  the  hole  another  major  fracture  zone 
was  encountered.  On  the  1st  of  June  flow  increased  from  80  to  240 
gpm  as  the  hole  was  advanced  from  1880  ft  to  its  1900  ft  final 
length.  Over  the  next  6  days  flow  decreased  slowly  to  225  gpm. 

A  series  of  flow  and  shut  in  experiments  was  conducted  on 
hole  Ul2t.04  UG-1  during  the  summer  and  fall  of  1973.  With  the 
flow  blocked  off,  pressure  in  the  hole  built  to  a  maximum  of  127 
psi  on  26  October.  This  corresponds  to  a  pressure  head  of  293  ft 
of  water.  The  hole  was  pressure  grouted  on  the  15th  of  November. 
Prior  to  grouting,  the  flow  was  generally  in  the  100  -  125  gpm 
range. 


The  flow  and  pressure  data  presented  in  Table  3.1 
indicate  a  range  of  peak  flows  in  the  various  exploration  holes 
of  from  100  to  240  gpm,  with  the  O12t.04  UG-1  drill  hole  producing 
the  greatest  flow.  Peak  shut  in  pressures  varied  from  42  to  165 
psi  corresponding  to  heads  of  between  97  and  381  ft.  This  large 
head  variation  supports  the  earlier  conclusions  of  Thordarson 
(1965)  that  the  elevation  of  fracture  saturation  varies 


32 


significantly  within  a  limited  area  due  to  the  generally  poor 
hydraulic  interconnectivity  of  the  fracture  zones. 

Indications  of  limited  hydraulic  interconnectivity  were 
provided  by  pressure  measurements  in  the  U12t.03  UG-1  hole  during 
drilling  of  the  Ul2t.05  UG-1  hole  and  by  pressure  measurements  in 
the  Ul2t.04  UG-1  hole  during  drilling  of  U12t.06  UG-1.  Upon 
completion  of  the  U12t.03  UG-1  hole  the  flow  was  blocked  off  and 
pressure  monitored  during  drilling  of  the  UG12t.05  UG-1  hole. 

The  latter  crossed  but  did  not  intersect  the  U12t.03  UG-1  hole  at 
a  length  of  1020  ft  in  U12t.05  UG-1.  As  shown  in  Table  3.2  (from 
Hoover,  1975),  as  the  length  of  U12t.05  UG-1  advanced  from  790 
to  1070  ft,  the  pressure  in  the  U12t.03  UG-1  hole  decreased  from 
160  to  about  130  psi.  Flow  out  of  the  Ul2t.05  UG-1  hole 
increased  from  1  gpm  to  about  100  gpm.  Hoover  concluded  that 
hydraulic  communication  between  the  two  holes  was  established 
along  two  faults  which  intersected  each  hole  near  the  crossing 
point. 


Following  completion  of  U12t.05  UG-1,  a  packer  was 
inserted  into  U12t.03  UG-1  just  beyond  the  intersection  of  the 
two  holes  to  prevent  the  large  flows  of  water  from  faults  and 
fractures  at  the  1800  to  2360  ft  depth  in  Ul2t.03  UG-1  from 
reaching  the  intersection  with  Ul2t.05  UG-1.  Total  flow  of  water 
out  of  Ul2t.05  UG-1  promptly  diminished  from  150  to  about  10  gpm, 
indicating  that  the  bulk  of  the  flow  from  U12t.05  UG-1  was 
actually  produced  by  the  water  bearing  fractures  intersected  by 
Ul2t.03  UG-1. 

Pressure  was  also  monitored  in  the  Ul2t.04  UG-1  hole 
during  drilling  of  the  adjacent  U12t.06  UG-1  hole.  Hoover  (1975) 
reports  that  there  was  some  correlation  between  activity  in  the 
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06  hole  and  pressure  in  the  04  hole.  A  stabilized  pressure  of 
120-122  psi  was  reached  in  the  04  hole  during  drilling  of  06. 
During  a  four  day  period  the  pressure  in  Ul2t.04  UG-1  dropped 
to  106  psi  while  hole  06  was  allowed  to  drain  freely.  Hole  06 
had  been  drilled  to  a  length  of  1400  ft  and  the  flow  rate  from 
hole  06  dropped  from  142  gpm  to  123  gpro  during  this  interval. 

The  pressure  in  hole  04  took  two  days  to  restabilize  at  120  -  122 
psi  once  drilling  was  resumed  in  hole  06. 

In  a  1974  memorandum.  Hoover  (1974)  predicted  initial 
inflows  of  water  into  the  Ul2t.03  HUSKY  PUP  main  drift  based  on 
flow  measurements  made  in  the  U12t.03  UG-3  exploration  hole.  The 
U12t.03  drift  followed  essentially  the  same  path  as  the  U12t.03 
UG-3  drill  hole.  Hoover  predicted  initial  flows  of  approximately 
60  and  125  gpm  through  fracture  zones  at  locations  in  the  Ul2t.03 
drift  corresponding  to  the  860  and  965  ft  ranges  in  the  U12t.03 
UG-3  drill  hole.  In  his  1975  memorandum  he  notes  that  neither  of 
these  anticipated  fracture  zones  were  encountered  in  the  Ul2t.03 
drift.  No  explanation  is  offered  for  this  discrepancy. 

The  most  significant  inflow  into  the  Ul2t.03  HUSKY  PUP 
drift  occurred  at  the  terminus  of  the  drift,  about  65  ft  short  of 
a  fracture  zone  which  was  predicted  to  flow  at  310  gpm  had  the 
drift  been  extended  that  far.  Hoover  believes  that  the  60  gpm 
flow  which  occurred  at  the  end  of  the  drift  came  from 
interconnecting  fractures  into  this  zone.  He  also  notes  that 
inflow  from  relatively  short  rockbolt  holes  attests  to  the  poor 
hydraulic  interconnectivity  in  this  area.  In  a  final 
observation.  Hoover  mentions  that  there  was  a  30  -  35  gpm  flow 
into  the  HUSKY  PUP  Bypass  Drift  extension  from  a  fracture  zone 
which  was  dry  in  an  adjacent  drift  only  10  ft  away.  This 
fracture  zone  had  produced  a  similar  initial  flow  rate  when 
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penetrated  by  the  adjacent  drift.  This  is  further  testimony  to 
zones  of  very  poor  hydraulic  interconnectivity. 

Hoover  (1975)  makes  the  following  conclusions  with  regard 
to  water  in  the  T  tunnel  complex: 

•  Most  of  the  water  flowing  from  the  exploratory  drill 
holes  is  evidently  from  poorly  connected  reservoirs 
in  fault  and  fracture  zones; 

•  The  pressure  response  in  the  U12t.03  UG-1  and  Ul2t.04 
(JG-1  drill  holes  to  activities  in  the  Ul2t.05  UG-1  and 
U12t.06  UG-1  drill  holes  indicates  that  some  individual 
fracture  or  fault  zones  are  interconnected  and  open  to 
water  flow  over  distances  of  at  least  several  hundred 
feet; 

•  Maximum  pressures  in  the  drill  holes  indicate  that  the 
pressure  head  in  the  fractures  varies  significantly 
within  this  relatively  small  area.  Maximum  pressures 
indicate  fracture  saturation  from  slightly  below  to  well 
below  the  top  of  the  zeolitized  tunnel  bed  tuffs. 
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Table  3.1.  Peak  flows  and  pressures  In  the  U12t  exploration  holes 
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Tcdale  3.2.  Pressure  response  in  Ul2t.03  UG-1  drill  hole 
to  drilling  in  U12t.05  UG-1  drill  hole  (from 
Hoover,  1975). 


U12t.03 

U12t.05 

UG-l 

UG-1 

pressure 

depth 

(psl) 

(feet) 

Remarks 

Maximum  pressure  U12t.05  UG-I 

not  drilled. 

Table  3.2  continued.  Pressure  response  in  Ul2t.03  UG-1 

drill  hole  to  drilling  in  Ul2t.05  UG-1 
drill  hole. 


Table  3.2  continxied.  Pressure  response  in  Ul2t.03  UG-1  drill 

hole  to  drilling  in  Ul2t.05  UG-1  drill 
hole . 


Dac« 

U12C.03 

UG-1 

pressure 

(pst) 

U12t.05 

UG-1 

depth 

(feet) 

U12C.05 

UC-1 

waterf low 
(Rpm) 

Remarks 

1973 

Aug.  27 

12  7 

2,352 

220 

Open  tools  at  200  ft. 

27 

120 

2,352 

195 

Open  tools  at  1,100  ft. 

28 

118 

2,352 

240 

Open  hole. 

28 

120 

2,352 

90 

All  tools  in  at  240  ft. 

30- 

Sept.  4 

Set  packer  In  U12t.03  UG-1 

at  1,026  ft. 

4 

142 

2,352 

54 

Open  tools  at  240  ft.  Swing  shift. 

5 

142-145 

2,352 

42 

Open  tools  at  240  ft.  Graveyard  shift. 

5 

142 

2,352 

3 

Open  tools  at  240  ft.  Day  shift. 

5 

147 

2,352 

3 

Open  tools  at  240  ft.  Swing  shift. 

6 

147 

2,352 

10 

Open  tools  at  240  ft.  Graveyard  shift. 

6 

143 

2,352 

12 

Open  tools  at  240  ft.  Day  shift. 

6 

150 

2,352 

10 

Open  tools  at  240  ft.  Swing  shift. 

6 

Packer 

in  U12t.03  UG-1  blew  up  at  2130  hrs. 

Sept.  7 

U12t.03  UG-1  open  or  drilling  on  packer.  Flou  in  U12t.05  UG-1  I 

to 

. 

Nov.  12 

used 

for  drilling. 

9 

U12t.03 

UC-1  shut— in  after 

drilling  out  packer.  75  psi 

U12t. 

03  UC-1  shut  in. 

12 

90 

2,352 

Ul2t.05  shut-in  at  120  psi. 

13 

90 

2,352 

Ul2t.05  shut-in  at  125  psi. 

Both  drill  holes  grouced  November  13,  1973 
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Figure  3.1.  Locations  of  major  txinnel  complexes.  Rainier  Mesa 
(courtesy  of  Dean  R.  Townsend,  Fenix  and  Scisson, 
Inc. ) . 


40 


iO 


41 


Figure  3.2.  Fault  locations  and  occurrence  of  fracture  water  in  E  Tunnel,  Rainier 
Mesa  (from  Thordarson,  1965). 
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Figure  3.3.  Tunneling  progress  and  total  flow  from  LOGAN  and  BLANCA  Drifts 
(from  Thordarson,  1965). 
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Figure  3.4.  Location  of  Ul2n.03  Drift  (from  Ege  et  al. ,  1980). 
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Ul2n.03  Drift 
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Figure  3.5.  Geology  of  the  Ul2n.03  drift  in  the  vicinity  of  the  Aqueduct  Syncline 
(from  Ege  et  al. ,  1980)  . 


Locations  of  T  Tunnel  exploratory  drill  holes. 


Hole  Length  (ft)  Total  Flow  (gpm) 


SECTION  4 


INFLOW  OF  WATER  FOLLOWING  UNDERGROUND  TESTS 


INTRODUCTION 


Inflow  of  water  following  underground  nuclear  detonations 
in  the  tunnel  bed  tuffs  of  Rainier  Mesa  has  not  been  a  major 
problem.  However*  following  the  MIGHTY  EPIC  event  in  the  Ul2n 
tunnel  complex,  there  was  flooding  of  access  drifts  and  severe 
flooding  of  .hardened  test  structures  in  the  vicinity  of  the 
detonation . 

There  are  a  number  of  possible  explanations  as  to  why 
such  flooding  has  not  been  a  more  common  occurrence.  Most 
obvious  is  the  location  of  many  of  the  underground  tests  in  areas 
where  the  level  of  fracture  saturation,  i.e.  the  easily  mobilized 
ground  water,  is  below  the  level  of  the  shot.  In  general,  tests 
in  the  G  and  B  tunnel  complexes  fall  into  this  category. 

A  second  explanation  for  lack  of  flooding  is  that  the 
level  of  fracture  saturation  varies  significantly  within  the 
tunnel  complexes  which  generally  lie  beneath  the  water  table.  As 
described  in  previous  sections,  there  are  apparently  "dry"  areas 
within  these  complexes  which  are  either  within  a  zone  of 
depressed  fracture  saturation  or  are  in  an  area  where  the  fractures 
drain  to  the  underlying  unsaturated  and  more  permeable  rock. 

Tests  within  such  areas  have  not  flooded  because  there  is  little  or 
no  fracture  water  to  mobilize. 
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A  third  explanation  £or  the  low  incidence  of 
flooding  is  the  possibility  that  the  underground  nuclear 
detonations  may  actually  induce  drainage  of  saturated  fractures 
into  the  underlying  unsaturated  rock.  Extensive  zones  of  block 
motion  occur  to  substantial  ranges,  in  some  instances  to  peak 
dynamic  stress  levels  as  low  as  0.1  to  0.25  Kbar,  surrounding  the 
underground  events.  Block  motions  are  permanent  displacements 
along  pre-existing  planes  of  weakness  such  as  faults  and  bedding 
planes.  These  have  been  well  documented  on  several  events,  e.g. 
Short  and  Kennedy  (1982),  Blouin  (1980).  Such  motions  surrounding 
and  beneath  a  shot  may  provide  drainage  paths  to  the  underlying 
more  permeable  rock,  thus  tending  to  dewater  rather  than  flood  a 
particular  site.  A  case  in  point  may  be  the  MIGHTY  EPIC  event  in 
T  Tunnel  adjacent  to  the  HUSKY  PUP  event.  Prior  to  the  test  many 
of  the  faults  and  bedding  planes  were  wet  or  seeping  water. 
Following  the  test  there  was  little  or  no  water  evident.  In 
addition,  the  presence  of  fracture  water  may  enhance  the 
occurrence  and  extent  of  block  motion  by  lubrication  of  fracture 
surfaces  and/or  by  preventing  mobilization  of  frictional 
resistance  during  dynamic  loading. 

A  final  factor  contributing  to  the  low  flooding  incidence  . 
is  the  fact  that  most  of  the  underground  nuclear  tests  are 
detonated  in  close  proximity  to  previous  tests  which  have  often 
already  loaded  the  surrounding  rock  to  very  significant  stress 
levels  and  may  have  already  mobilized  fracture  water  in  the  area 
and/or  drained  the  site. 

For  whatever  reasons,  flooding  on  the  Rainier  events  has 
not  been  widespread.-  The  flooding  on  MIGHTY  EPIC,  however,  stands 
as  a  warning  that  the  problem  cannot  be  ignored.  In  fact,  in  a 
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different  geologic  setting  where  rock  permeability  is  higher  than 
that  in  the  tuffs  and  where  the  water  table  is  not  perched,  as  at 
Rainier  Mesa,  flooding  similar  to  that  experienced  on  MIGHTY  EPIC 
might  be  the  rule  rather  than  the  exception. 

Two  mechanisms  are  often  put  forth  to  explain  the 
generation  of  water  from  underground  explosions  in  the  NTS  tuffs. 
In  the  first,  water  stored  in  fracture  zones  above  and  adjacent 
to  the  underground  works  is  released  by  explosive  disturbance  of 
previously  impermeable  fracture  zones  leading  to  the  tunnels. 

Such  disturbances  may  be  in  the  form  of  block  motion  or  rebound, 
whereby  there  is  dilatency  in  the  fractures  allowing  the 
water  to  flow  rapidly  into  the  tunnels. 

The  second  mechanism  is  akin  to  the  consolidation  process 
in  soil  mechanics.  In  this  process  a  residual  stress  field 
resulting  from  the  explosive  detonations  is  imposed  on  the 
saturated  rock.  The  residual  stresses  squeeze  both  the. 
interstitial  water  and  fracture  water  from  the  surrounding  rock 
into  the  underground  works.  In  this  process  the  amount  of  water 
produced  and  the  rate  of  flow  are  dependent  on  the  permeability 
of  both  the  interstitial  rock  and  fractures,  the  magnitude  of  the 
residual  stresses,  and  the  net  porosity  and  mechanical 
properties  of  the  rock. 

There  is  solid  experimental  evidence  for  both 
of  the  above  water  production  mechanisms  on  underground  nuclear 
and  HE  tests  beneath  Rainier  Mesa.  Following  reentry  to  MIGHTY 
EPIC  extensive  flooding  was  observed,  with  water  still  flowing 
from  several  faults  which  had  been  dry  prior  to  the  test  and 
other  evidence  of  flow  from  faults  which  had  been  displaced  by 
the  shot.  On  the  HE  shot  ONETON,  conducted  in  G  Tunnel, 
generally  considered  to  be  above  the  level  of  fracture 
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saturation,  a  substantial  amount  of  water  was  produced  from  the 
walls  and  floor  of  an  open  drift  located  at  a  range  where  the 
peak  dynamic  stress  was  only  0.12  Kbar.  This  appeared  to  be 
interstitial  water  squeezed  from  the  intact  rock  by  the  elevated 
residual  stresses  generated  by  the  explosion.  Water  flowed  for 
hundreds  of  hours  following  the  detonation  and  was  accompanied  oy 
continually  decreasing  residual  stresses  within  the  rock  mass 
adjacent  to  the  drift. 

The  following  subsections  describe  available  details  on 
the  MIGHTY  EPIC  and  ONETON  water  production. 


MIGHTY  EPIC 

MIGHTY  EPIC  was  a  low  yield  (i.e.  less  than  20  KT) 
contained  nuclear  event  fielded  in  the  Ul2n.l0  drift  beneath 
Rainier  Mesa.  The  location  of  MIGHTY  EPIC  is  shown  on  the  plan 
view  of  Figure  4.1.  The  working  point  was  located  in  Tunnel  Bed 
3  whi.le  most  of  the  test  structures  and  the  line  of  site  pipe  were 
primarily  in  the  zeolitized  tuffs  of  Tunnel  Bed  4.  According  to 
Townsend  (1984)  there  was  virtually  no  water  inflow  during  mining 
of  the  MIGHTY  EPIC  line  of  site  and  bypass  drifts.  The  only 
signs  of  water  were  a  few  damp  fracture  surfaces.  Also,  unlike 
the  majority  of  tests,  the  MIGHTY  EPIC  works  were  in  virgin 
ground  which  had  not  been  subjected  to  high  previous  stresses  from 
adjacent  shots. 

A  plan  view  of  the  MIGHTY  EPIC  event  is  shown  in  Figure 
4.2.  The  working  point  is  at  the  west  end  of  the  U12n.l0  main 
drift.  This  drift  houses  the  line  of  sito  (LOS)  pipe  which 
contained  targets  for  exposure  to  radiation  from  the  MIGHTY  EPIC 
device.  The  LOS  pipe  is  simply  a  large  diameter  tapered  steel 
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tube  with  mounts  for  the  targets  located  at  prescribed  ranges. 
These  targets  are  shielded  from  blast  and  debris  by  two  large 
gas  driven  sliding  doors  at  locations  labeled  DAC  1  and  DAC  2. 
Beyond  the  DACs  is  a  third  closure,  called  the  Tunnel  and  Pipe 
Seal,  or  TAPS,  which  was  located  approximately  100  ft  east  of  the 
DAC  2.  The  radiation  targets  are  located  to  the  east  of  the  TAPS 
at  various  ranges  along  the  LOS  pipe. 

The  Bypass  Drift  runs  parallel  to  the  LOS  Drift  and 
provides  access  to  the  working  point  room  and  Interface  Drift 
during  construction  of  the  LOS  pipe.  Prior  to  the  test  the 
Bypass  Drift  is  grouted  closed  to  the  range  indicated  by  the  end 
of  stemming  between  the  B  and  C  Structures  Drifts.  Beyond  the 
end  of  stemming  the  Bypass  drift  is  open  and  is  reinforced  with  a 
standard  rockbolt,  wire  mesh  and  gunite  lining. 

Running  SW  off  the  Bypass  Drift  are  th<a  A,  B,  and  C 
Structures  Drifts.  These  contained  horizontal  cylindrical 
hardened  structures  at  ranges  of  290,  400  and  600  ft  from  the 
working  point.  Various  structural  concepts,  sizes  and  strengths 
were  tested  in  these  drifts. 

The  Interface  Drift  running  north  from  the  working  point 
provided  access  to  instrumentation  holes  drilled  vertically 
downward  through  the  tunnel  bed  tuff  to  the  underlying  quartzite. 
Instrumentation  in  these  holes  was  designed  to  measure  block 
motion  betwee"  the  tuff  and  underlying  quartzite.  The  Interface 
Drift  was  also  grouted  shut  prior  to  the  test. 

The  post-shot  observations  of  water  on  the  MIGHTY  EPIC 
event  are  based  on  two  memos  from  Dean  R.  Townsend  (1976,  1977) 
and  on  informal  interviews  with  him  (1984)  in  which  he  supplied 
additional  details  from  memory.  He  reentered  the  MIGHTY  EPIC 
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Bypass  Drift  approximately  2  days  after  the  event  through  the 
overburden  plug.  This  is  an  additional  safety  plug  to  keep  any 
radiation  which  might  escape  near  the  source  region  from 
entering  the  tunnel  complex.  It  is  located  just  beyond  the  end 
of  the  MIGHTY  EPIC  Bypass  Drift  at  a  range  of  about  2000  ft  from 
the  working  point.  On  reentry  the  Bypass  Drift  was  flooded 
from  the  end  of  stemming  all  the  way  to  the  overburden  plug/  a 
distance  of  approximately  1600  ft.  There  was  approximately  6 
inches  of  water  on  the  floor  of  the  drift,  the  water  being 
somewhat  deeper  near  the  overburden  plug  and  shallower  near  the 
end  of  stemming  due  to  the  gradient  of  the  tunnel.  Taking  the 
average  tunnel  width  as  11  ft,  it  is  estimated  that  there  was 
probably  in  excess  of  50,000  gal  of  water  in  the  Bypass  Drift,  A 
similar  amount  was  probably  contained  in  the  LOS  Drift. 

Water  flowing  from  two  faults  appears  to  have  been  the 
primary  source  of  this  flooding.  A  small  reverse  fault  in  the 
Bypass  Drift  at  a  range  of  870  ft  from  the  working  point  was 
flowing  at  3  -  5  gpm  several  days  after  the  event.  A  larger 
normal  fault  in  the  LOS  Drift  at  a  range  of  940  ft  from  the 
working  point  was  emitting  5  to  10  gpm  at  this  time.  There  was 
no  evidence  of  significant  block  motion  on  these  faults,  though 
the  gunite  lining  had  been  spalled  from  the  walls,  evidently  by 
the  flowing  water.  These  faults  continued  to  flow  for 
approximately  20  days  after  reentry.  Since  there  was  no  flow 
from  these  faults  during  mining,  the  MIGHTY  EPIC  detonation  must 
have  either  opened  drainage  to  these  faults  (if  they  were 
already  open),  or  opened  these  faults  sufficiently  to  allow 
drainage  from  saturated  zones  within  them  and/or  from  within 
other  saturated  fractures  connected  to  them. 

During  the  initial  reentry,  water  was  observed  in  the  LOS 
pipe  between  the  DAC  1  closure  and  the  TAPS  and  water  was 
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observed  leaking  from  around  DAC  1.  The  source  of  this  water  was 
believed  to  be  inflow  from  Fault  5,  which  ran  NS  and  intersected 
the  LOS  Pipe  between  the  two  DACs.  Block  motion  occurred  along 
Fault  5,  with  about  1.5  ft  of  lateral  motion  in  the  area  of  its 
intersection  with  the  Bypass  and  LOS  Drifts.  Figure  4.3  shows 
the  location  of  major  faults  and  bedding  planes  within  the  MIGHTY 
EPIC  structures  region  with  the  direction  and  magnitudes  of  block 
motion  indicated  on  the  planes  of  weakness  which  were  activated 
by  the  detonation.  Townsend  (1977)  reports  that  post-test  mining 
eventually  provided  access  to  the  LOS  Pipe  between  the  two  DACs. 
The  LOS  Pipe  was  partially  crushed  and  torn  open  at  its 
intersection  with  Fault  5  adjacent  to  the  DAC  2  closure.  A  view 
of  the  pipe  at  this  location  is  shown  in  Figure  4.4.  The 
interior  of  the  LOS  pipe  between  the  two  DACs  appeared  to  have 
been  full  of  water  shortly  after  the  test,  as  evidenced  by  mud 
and  water  marks  within  the  pipe.  This  water  slowly  leaked  out 
past  the  DAC  1  over  a  period  of  days  following  the  test. 

Several  centimeters  of  pink  mud  were  deposited  in  the  LOS 
pipe  between  the  DACs.  This  mud  consisted  of  finely  powdered 
zeolitized  tuff  which  appeared  to  be  identical  to  a  1  cm  wide 
seam  of  pale  pink  fault  gouge  lining  Fault  5  in  this  region. 

Based  on  this  evidence  Townsend  (1977)  concludes  that  movement 
along  Fault  5  provided  a  channel  for  post-test  flow  into  the  LOS 
Pipe.  The  actual  source  of  this  water  was  probably  fracture 
water  trapped  in  Fault  5  or  in  other  fractures  connected  to  Fault 
5. 


Post-test  mining  and  reentry  to  the  hardened  experimental 
structures  in  drifts  A  and  B  revealed  that  these  structures  were 
full  of  water.  As  shown  in  Figure  4.3,  the  block  motion  along 
Fault  5  also  intersected  the  B  Structures  Drift  at  which  point 
there  was  a  relative  lateral  displacement  of  1.3  ft.  This  motion 
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heavily  damaged  the  structure  at  the  intersection,  buckling  and 
tearing  the  steel  liner  and  evidently  allowing  water  from  Fault  5 
to  completely  flood  the  structure.  Views  of  the  damage  from  the 
Fault  5  relative  displacement  are  shown  in  Figure  4.5, 

Reentry  to  the  A  Structures  Drift  revealed  that  one  of 
the  structures  had  been  torn  open  in  the  vicinity  of  Fault  7, 
thus  providing  an  entry  point  for  the  water  which  flooded  that 
drift.  No  significant  block  motion  was  indicated  on  Fault  7; 
however,  so  the  path  and  source  of  the  water  filling  A  Drift  is 
uncertain. 

In  summary,  many  tens  of  thousands  of  gallons  of  water 
flooded  the  MIGHTY  EPIC  works  at  a  number  of  locations  following 
the  detonation.  The  sources  of  most  of  this  water  were  probably 
saturated  faults  and  fracture  zones  which  were  tapped  or 
mobilized  by  the  stresses  and  motions  generated  by  the  explosion. 
In  several  locations  flow  evidently  occurred  along  faults  which 
were  displaced  more  than  a  foot  by  the  test.  Steel  lined  hardened 
structures  were  completely  flooded  by  water  entering  through 
tears  in  the  liners.  Significant  flows  of  water  were  observed 
from  faults  at  great  ranges  (nearly  1000  ft)  from  the  working 
point.  Stress  levels  and  motions  at  these  ranges  are  very  low, 
well  below  levels  normally  associated  with  damage  to  even 
minimally  lined  tunnels.  The  MIGHTY  EPIC  drifts  were  dry  during 
and  after  mining;  the  substantial  flooding  which  occurred 
resulted  solely  from  disturbance  by  the  explosion. 


ONETON 


The  ONETON  HE  event  in  G  Tunnel  produced  water  in  a 
manner  which  contrasts  with  that  described  previously  on  MIGHTY 
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EPIC.  Instead  of  mobilizing  fracture  water,  the  ONETON  event 
forced  interstitial  water  out  of  the  voids  in  the  saturated  tuff. 
The  residual  stresses  and  stress  gradients  produced  by  the 
detonation  apparently  generated  pore  pressures  within  the  rock 
which  caused  the  pore  water  to  migrate  away  from  the  higher 
stressed  regions  near  the  shot  point  and  to  flow  into  an  open 
drift  some  distance  away. 

The  layout  of  the  ONETON  event  is  shown  in  the  plan  view 
of  Figure  4.6.  A  2000  lb  TNT  sphere  was  detonated  at  the  end  of 
a  fully  stemmed  dogleg  drift.  Stress  and  motion  gages  to 
actively  monitor  total  dynamic  and  residual  stresses  and  ground 
motions  were  placed  at  locations  1  through  12.  Instrumentation 
cables  were  fed  into  the  open  drift  adjacent  to  the  bulkhead  at 
the  end  of  stemming.  Smith  (1984)  described  the  ONETON  setting 
in  G  Tunnel  as  saturated  (about  98%)  and  free  of  faults  and 
joints.  No  water  was  encountered  during  the  ONETON  mining. 

Peak  dynamic  radial  stress  as  a  function  of  range  is 
plotted  in  Figure  4.7.  The  ONETON  data  are  a  good  match  to  the 
scaled  peak  stress  attenuation  from  a  previous  series  of  64  lb  HE 
shots  in  this  material.  The  ONETON  peak  stress  data  are  also  in 
reasonably  good  agreement  with  scaled  peak  stress  data  from  the 
contained  nuclear  events. 

Smith's  (1983)  plot  of  residual  radial  stress  as  a 
function  of  range  for  6  HE  shots  in  nearly  saturated  tuff  is 
shown  in  Figure  4.8.  The  PUFF  TOO  1000  lb  data  and  the  ONETON 
2000  lb  data  have  been  scaled  by  the  cube  root  of  their  yields  to 
be  consistent  with  the  64  lb  data.  There  is  considerable  scatter 
in  the  residual  stress  data,  though  deletion  of  the  two  low  data 
points  from  the  RS  14  event,  which  Smith  notes  are  in  an 
unusually  soft  tuff  layer,  considerably  improves  the  picture. 
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Note  that  residual  stresses  are  considerably  less  than  the  peak 
stresses  at  any  given  range  and  that  they  tend  to  decay  less 
rapidly  than  the  peak  stresses. 

The  residual  stress  field  primarily  results  from 
formation  of  a  permanent  cavity  by  the  detonation.  Close  to  the 
explosion  the  high  amplitude  dynamic  stresses  cause  severe 
plastic  deformation  of  the  rock  and  a  strong  outward  thrust  of 
material  as  the  cavity  around  the  detonation  expands  dynamically. 
Following  the  maximum  dynamic  expansion  of  the  cavity  there  is  a 
rebound  caused  by  the  high  imbalanced  stresses  in  the  surrounding 
material.  The  cavity  undergoes  a  partial  compression,  finally 
coming  into  equilibrium  with  stresses  in  the  surrounding  rock. 

The  resulting  equilibrium  stresses  are  very  high  in  the  severely 
deformed  rock  surrounding  the  cavity,  dropping  monotonically  to 
lower  values  in  the  less  severely  distorted  rock  further  from  the 
cavity.  Smith  notes  that  the  ONETON  cavity,  measured  after 
reentry,  was  quite  uniform,  with  an  average  radius  of  4.4  ft. 

The  decay  of  residual  stress  with  time  on  the  ONETON 
event  is  shown  in  Figure  4.9  for  a  period  of  about  17  hours  after 
the  event.  Residual  stresses  decayed  to  approximately  half 
their  immediate  post-test  values  over  this  time  span.  Similar 
data  were  recorded  on  the  64  lb  events.  Smith  suggests  that  the 
decay  of  residual  stress  with  time  is  related  to  migration  of  pore 
water  induced  by  the  residual  stress  gradient.  He  sites  seepage 
of  water  into  the  open  drift  pictured  in  Figure  4.6  as  evidence 
of  this  migration. 

Following  the  ONETON  detonation,  water  was  observed 
seeping  out  of  the  face,  ceiling,  walls,  and  floor  of  the  open 
drift  (Smith,  1984).  The  face  was  located  at  a  range  of  55  ft 
from  the  working  point.  Following  the  detonation  damp  zones 
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appeared  on  the  walls  and  ceiling.  These  became  increasingly 
moist  until  finally  water  began  running  off  the  ceiling  and  down 
the  walls  and  collected  on  the  floor  of  the  drift.  During  the 
first  121  hours  after  the  shot,  215  gallons  of  water  seeped  out. 
During  the  following  116  hours,  an  additional  335  gallons  seeped 
from  the  face  and  sides  of  the  drift.  Smith  notes  that  during 
the  reentry  mining,  no  distinct  fracture  zones  or  other  water 
bearing  aquifers  were  encountered,  suggesting  that  pore  water  was 
migrating  away  from  the  high  residual  stress  gradients.  As  the 
pore  water  migrates  out  of  the  higher  stressed  rock,  pore 
compression  or  partial  collapse  will  occur  accompanied  by  a 
relaxation  of  both  pore  pressures  and  total  stresses  in  this 
region.  This  mechanism,  which  is  similar  to  the  consolidation 
process  in  saturated  fine  grained  soils,  would  account  for  the 
decaying  residual  stresses  recorded  in  Figure  4.9.  A  preliminary 
two  phase  calculation  of  the  ONETON  experiment  is  described  in 
Section  7. 
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Figure  4.1, 


Location  of  MIGHTY  EPIC  event 
(from  Short  and  Kennedy,  1982) 
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Figure  4.3.  Location  arid 
motion  (from 


C  Drift 


U 


tude  of  MIGHTY  EPIC  block 
;  and  Kennedy,  1982) . 


6] 


UNCLASSIFIED 


a)  Water  and  damage  due  to  movement 
along  fault  '5  looking  southwest. 


h)  Damage  due  to  movement  along 
*ault  -5  U'or.ino  noi't'ieas I . 


Figtire  4.5.  Damage  to  B  Structures  Drift  from  relative 
displacement  along  Fault  5  (from  Short  and 
Kennedy,  1982). 
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Figure  4.6.  ONETON  test  layout  -  plan  view  (from  Smith,  1983). 


STRESS  *’  kbar 


Figure  4.7.  ONETON  peak  dynamic  radial  stresses  (from 
Smith,  1983). 
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Figure  4.9.  Decay  of  ONETON  residual  stresses  (from  Smith,  1984). 


SECTION  5 


CONSL  CODE  DESCRIPTION 


INTRODUCTION 

The  numerical  parametric  calculations  in  this  initial 
study  were  conducted  using  the  one  dimensional/  axisymmetric,  two 
or  three  phase,  quasi-static  finite  element  code  CONSL.  CONSL  is 
an  adaption  of  the  general  two  dimensional  quasi-static  analysis 
program  QSAP  written  and  described  by  Kim  (1982).  CONSL  can  be 
used  to  solve  both  uniaxial  strain  and  axisymmetric  flow  and 
consolidation  problems.  It  was  written  in  standard  FORTRAN  and 
can  even  be  run  on  microcomputers  using  the  CP/M  operating 
system.  CONSL  can  model  nonlinear  material  behavior  and 
calculate  flow  rates  and  consolidation  in  multilayered  media. 
Imposed  loading  conditions  can  include  specified  total  stresses, 
pore  pressures  and  skeleton  displacements.  CONSL  can  also  solve 
three  phase  problems  in  which  the  degree  of  saturation  is  over 
85%.  CONSL  was  checked  against  Terzaghi's  closed  form  solution 
for  one  dimensional  consolidation.  There  was  excellent  agreement 
between  the  closed  form  solution  and  the  CONSL  results.  A  copy 
of  CONSL  is  included  in  Appendix  D.  This  version  is  written  in 
FORTRAN  77  for  an  HPIOOO  system. 

In  the  next  subsection  a  description  of  the  input  to  the 
CONSL  program  is  given  which  serves  to  briefly  describe  the 
features  of  the  program  as  well  as  to  provide  a  guide  to  its 
use.  The  mathematical  finite  element  formulations  used  in  CONSL 
are  described  in  Appendix  C.  Appendix  A  presents  derivations  of 
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the  undrained  properties  of  saturated  media  in  plane  strain. 
These  were  used  to  check  the  initial  loading  conditions  computed 
using  CONSL.  Appendix  B  is  a  derivation  of  one  of  the  two  flow 
algorithms  used  in  CONSL. 

DESCRIPTION  OF  INPUT  DATA 


Input  must  be  in  terms  of  consistent  units  such  as  SI, 
English,  etc. 

Card  1 


Problem  Title  (80  characters) 


Card  2 

RI,  RO,  AP,  NUMEL,  NKF,  NEF,  NTF,  IPLANE,  NVF ,  NDF  (3E10.0,7I5) 

RI  =  radius  of  the  tunnel 
RO  =  radius  of  remote  boundary 
AP  =  mesh  growth  factor 


If 

AP 

= 

0.5,  element  sizes 

are  constant 

If 

AP 

< 

0.5,  element  sizes 

become  larger 

toward  the 

remote  boundary 

If 

AP 

> 

0.5,  element  sizes 

becomes  smaller  toward  the 

remote  boundary 

For 

AP 

= 

0.25,  half  of  the 

total  elements 

are  located 

within  the  first  25%  of  the  total  distance 
to  the  remote  boundary 
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NUMEL  =  total  number  of  elements 


NKF  = 


NEF  = 


NTF  = 


I  PLANE  = 


NVF  = 


NDF  = 


The  size  of  the  first  element  ai,  is  given  by 

=  RO  -  RI  r  1  NUMEL  -  1  , , 

1  NUMEL  NUMEL -  (5-1) 

0  linear  elastic  skeleton 
1  nonlinear  skeleton 

0  saturated  skeleton 
1  partially  saturated  skeleton 

0  constant  time  steps 
1  variable  time  steps 

0  1-0  cylindrical  symmetry 

1  1-D  plane  strain 

0  flow  volume  calculated  from  pore  pressure 
gradient  on  first  element 
1  flow  volume  based  on  Equation  B-19  as 
described  in  Appendix  B 

0  decoupled  two  phase  material  model  as 
described  in  Appendix  A 
1  fully  coupled  two  phase  material  model  as 
described  in  Appendices  A  and  C. 
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Card  3 


CM,  V,  POR,  PK,  (E10.0,2F10.0,E10.0) 

CM  =  constrained  skeleton  modulus 
V  =  Poisson's  ratio  of  skeleton 
POR  =  porosity 

PK  =  coefficient  of  permeability 

Card  £  (for  NEF  =  1) 

SO,  STAW  (2E10.0) 

SO  =  degree  of  saturation 

STAW  =  pressure  difference  between  pore  water  and  pore 
air  (see  Kim,  1982) 


Card  5 
TEND  (ElO.O) 

TEND  =  maximum  calculation  time 
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Card  6  (for  NTF  =  0) 

DT  (ElO.O) 

DT  =  constant  time  step 

Card  Group  6  (for  NTS  =  1) 

NDT(I5)  . 

NCL(l),  DTT(l)  I 

NCL(2),  DTT(2)  / 

.  >  NDT  Cards  with  (15,  ElO.O) 

NCL(NDT),  DTT(NDT)  ' 

NDT  =  number  of  different  time  steps 

NCL  =  number  of  cycles  having  given  time  step  DTT 

DTT  =  duration  of  the  given  zime  step 

Card  2. 

SRI,  STI,  PI  (3E10.0) 

SRI  =  initial  effective  radial  stress  (compression  is 
negative ) 

STI  =  initial  effective  tangential  stress 
PI  =  initial  pore  fluid  pressure 
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Card  8  (left  end  boundary  conditions,  see  Figure  7.1) 


JSDL,  JDCL,  TNFL,  BPPL,  SSDL  (215,  3E10.0) 

JSDL  =  0  free  skeleton  boundary 

=  1  fixed  or  specified  skeleton  boundary 

JDCL  =  0  permeable  boundary 

=  1  impermeable  boundary 

TNFL  =  applied  total  stress  on  boundary 

BPPL  =  applied  pore  pressure  at  boundary 

SSDL  =  specified  boundary  displacement  (for  JSDL  =  1) 

Card  9^  right  end  boundary  conditions) 

JSDR,  JDCR,  TNFR,  BPPR,  SSDR  (2I5,3E10.0) 

JSDR  =  0  free  skeleton  boundary 

=  1  fixed  or  specified  skeleton  boundary 

JDCR  =  0  permeable  boundary 

=  1  impermeable  boundary 

TNFR  =  applied  total  stress  on  boundary 
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BPPR  =  applied  pore  pressure  at  boundary 

SSDR  =  specified  boundary  displacement  (for  JSDL  =  1) 


SECTION  6 


DEVELOPMENT  OF  RESIDUAL  STRESSES 


INTRODUCTION 

Blast  loadings  will  cause  significant  flow  of  inter¬ 
stitial  pore  water  from  in  situ  rock  masses  only  if  the  dynamic 
stresses  produce  long  term  residual  stresses  and/or  pore 
pressures.  Significant  pore  water  flow  cannot  develop  during  the 
dynamic  portion  of  the  loading  because  the  dynamic  loading  is  of 
extremely  short  duration  with  respect  to  the  rock  permeability. 

Two  mechanisms  for  development  of  residual  stresses  and 
pore  pressures  are  described  in  this  section.  The  first  is  the 
primary  cause  of  residual  stresses  surrounding  underground 
explosions.  The  second  is  the  cause  of  residual  pore  pressures 
beneath  near  surface  explosions. 

It  is  well  known  that  the  underground  nuclear  shots  in 
Rainier  Mesa  tuff  produce  an  extensive  zone  of  residual  stresses 
around  the  cavities  formed  by  the  detonations.  These  have  been 
calculated  by  Patch  (1984),  Rimer  and  Friedman  (1978)  and  others. 
In  addition,  measurements  of  long  term  residual  stresses 
surrounding  a  nuclear  detonation  in  Rainier  Tuff  have  been 
reported  by  Ellis  and  Kibler  (1983).  While  these  do  not 
precisely  match  the  calculations,  similar  stress  patterns  were 
measured  and  development  of  the  residual  stress  field  was 
documented.  Finally,  the  contained  HE  tuff  experiments  described 
by  Smith  (1983)  produced  residual  stress  fields  similar  to  those 
formed  by  the  nuclear  detonations. 
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As  described  in  Section  4,  the  residual  stress  field 
surrounding  an  explosively  formed  underground  cavity  results  from 
the  severe  plastic  distortion  of  the  rock  in  the  immediate 
vicinity  of  the  cavity.  The  rock  surrounding  the  cavity 
undergoes  severe  plastic  deformations  as  it  is  thrust  outward  by 
the  force  of  the  explosion.  The  permanent  cavity  expansion  locks 
very  high  compressive  stresses  into  this  rock  and  produces 
residual  equilibrium  stresses  extending  many  radii  outward  from 
the  cavity. 

The  second  residual  stress/pore  pressure  formation 
mechanism  is  due  to  the  hysteretic  nature  of  the  saturated  rock 
skeleton.  Nonrecoverable  skeleton  strain  produced  by  the  dynamic 
loading  generates  excess  residual  pore  pressures.  In  essence/  a 
portion  of  the  in  situ  stress  carried  by  the  rock  skeleton  before 
the  dynamic  loading  is  transferred  to  the  pore  water  following 
the  hysteretic  unloading.  This  second  mechanism  would  be 
expectea  to  occur  following  detonation  of  a  nuclear  surface 
burst. 


RESIDUAL  STRESS  FROM  CAVITY  EXPANSION 

Essoglow  and  Rogich  (1965)  presented  a  simple  method  for 
computing  the  residual  displacements  and  strains  surrounding  a 
contained  detonation.  In  their  formulation  a  spherical  chamber 
surrounding  the  explosive  device  is  expanded  outward  to  form  a 
spherical  cavity  of  final  radius  r^,.  The  volume  of  material 
between  the  initial  spherical  shell  and  the  final  cavity  shell  is 
assumed  to  be  redistributed  throughout  the  surrounding  rock. 

This  redistribution  of  material  is  assumed  to  occur  with  no 
volume  reduction  (noncompressibility)  and  with  no  loss  of 
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material  from  vaporization.  In  reality,  both  the  pre-shot  cavity 
volume  and  vaporization  volume  are  negligible  with  respect  to  the 
final  cavity  volume  for  an  underground  nuclear  explosion.  Based 
on  the  above  assumptions,  the  permanent  displacement,  6  ,  at  range 
r  is  given  by 


6 


1/3 


r 


(6-1) 


For  ranges  of  two  cavity  radii  and  beyond.  Equation  6-1  is 
closely  approximated  by 


6  * 


(6-2) 


3r 


The  permanent  radial  strain  e-,  is  given  by 


^r  = 


773 


+  1 


(6-3) 


For  ranges  beyond  two  cavity  radii  this  can  be  approximated  by 


2r. 


^r  = 


3r' 


(6-4) 


Even  though  the  assumptions  governing  the  derivation  of 
the  above  equations  are  quite  simplified,  they  appear  to  be  in 
good  agreement  with  test  data.  Figure  6.1  is  a  plot  of  Equation 
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6-2  compared  to  the  permanent  displacement  data  on  the  Rainier 
event  presented  by  Diment  et  al.  (1959).  Using  material  models 
for  appropriate  rock  types,  it  is  possible  to  compute  residual 
stresses  from  Equations  6-3  or  6-4. 

More  sophisticated  residual  stress  calculations  have 
become  available  in  recent  years  based  on  finite  difference 
calculations  which  model  the  entire  ground  motion  and  stress  time 
response  surrounding  underground  tests.  Examples  of  such 
calculations  are  given  in  Rimer  and  Friedman  (1978)  and  Patch 
(1984).  Figure  6.2,  from  Patch,  shows  computed  radial  and 
tangential  residual  stresses  as  a  function  of  range  normalized  by 
the  cavity  radius  for  four  recent  underground  nuclear  events  in 
Rainier  Mesa  tuff.  Significant  radial  and  tangential  stresses 
are  induced  to  ranges  of  many  cavity  radii.  In  the  calculations 
of  Section  7  a  residual  stress  of  1000  psi  (6.9  MPa)  was  assumed 
which  corresponds  to  the  radial  residual  stress  at  between  5  and 
6  cavity  radii  for  the  calculations  shown  in  Figure  6.2. 

RESIDUAL  EXCESS  PORE  PRESSURE  FROM  SURFACE  BURST  LOADING 

Residual  excess  pore  pressures  can  be  induced  by 
nonrecover able  skeleton  deformation  in  saturated  materials. 

While  the  skeleton  tends  to  only  partially  recover  from  a  dynamic 
loading,  the  pore  water  is  essentially  elastic  and  tends  to 
rebound  fully.  The  result  is  that  some  of  the  initial  effective 
stress  originally  carried  by  the  skeleton  is  taken  by  the  pore 
water.  In  the  limit,  all  of  the  effective  skeleton  stress  is 
transferred  to  the  pore  water  and  a  state  of  liquefaction  is 
achieved  (see  Kim  and  Blouin,  1984). 
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In  this  subsection  the  residual  excess  pore  pressures 
generated  by  simple  one  dimensional  uniaxial  loadings  are 
derived.  In  the  case  of  a  nuclear  surface  burst  the  actual 
loadings  are  considerably  more  complex  than  the  one  dimensional 
analog  discussed  here. 

To  simplify  the  complicated  two  phase  response,  the 
following  assumptions  have  been  made: 

•  Motions  are  vertically  one  dimensional.  In  practice, 
uniaxial  one  dimensional  loadings  adequately  simulate 
the  response  of  the  underlying  geology  to  airblast 
loadings. 

•  The  skeleton  uniaxial  stress-strain  curve  is  approx¬ 
imated  by  two  linear  slopes  consisting  of  a  loading 
constrained  modulus  (Mg'^)  and  an  unloading 

constrained  modulus  (Mg^).  This  bilinear 
hysteretic  skeleton  model  is  shown  in  Figure  6.3. 

•  The  solid  grains  and  pore  water  are  linearly  elastic. 

•  Drainage  is  not  permitted  during  the  dynamic  loading 
and  unloading. 


Partially  Coupled  Model 

If  volume  change  in  the  soil-water  mixture  due  to 
effective  stress  on  the  individual  grains  is  neglected,  the 
undrained  constrained  modulus  of  the  bulk  mixture,  Mp,  is 
given  by 
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M  =  M  - 
P  s 


K  K 
m  s 

K 


+  K 


rti 


(6-5) 


where  Mg  =  constrained  modulus  of  solid  skeleton 
K_  =  bulk  modulus  of  solid  skeleton 


K  K 


= 

m 


"a  " 

"w  = 
n  = 


K 


K  )  '  "'ixture  modulus 

g 


w 


bulk  modulus  of  solid  grains 
bulk  moduxus  of  pore  water 
porosity 


The  undrained  partially  coupled  modulus  was  derived  by  Blouin  and 
Kim  (1984)  and  is  included  in  Table  A.l.  Due  to  the  hysteretic 
nature  of  the  solid  skeleton,  the  unloading  modulus  is  greater 
than  the  loading  modulus.  The  loading  undrained  constrained 
modulus  of  the  bulk  mixture  is  given  by 


"p"'  =  '^s'"  -  ^  ^m 

g 

where  Mg^  and  Kg^  are  the  loading  constrained  and 

bulk  moduli  of  the  solid  skeleton  respectively.  The  unloading 

undrained  constrained  modulus  of  the  bulk  mixture  is  given  by 


M 

P 


U 


K  K 
m  s 


U 


K 

g 


+  K 


m 


(6-7) 
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where  Mg^  and  Kg^  are  the  unloading  constrained 
and  bulk  moduli  of  the  solid  skeleton  respectively. 


Figure  6.4  shows  volumetric  strain  change  during  the 
dynamic  loading  and  unloading.  Volumetric  strain  at  the  end  of 
the  loading  is  calculated  by 


(6-8) 


where  is  the  dynamic  total  vertical  stress  increment. 

The  recovered  volumetric  strain  at  the  end  of  the  unloading  is 
calculated  by 


(6-9) 


The  pore  pressure  increment  at  the  end  of  the  loading  is  related 
to  the  volumetric  strain  increment  by 


ATTp,^  =  K  Ae  ^ 
D  m  v 


(6-10) 


and  the  pore  pressure  drop  at  the  end  of  the  unloading  is  related 
to  the  recovered  volumetric  strain  by 


A  U  V  A  U 


(6-11) 


80 


The  residual  excess  pore  pressure  is  defined  as 


Substitution  of  Equations  6-10  and  6-11  into  Equation  6-12  gives 

-  AC^°)  (6-13) 

Using  Equations  6-8  and  6-9,  the  volumetric  terms  in  Equation 
6-13  can  be  expressed  in  terms  of  stresses  and  moduli. 

(6-14) 

”p  / 

Or,  the  ratio  of  excess  pore  pressure  to  the  dynamic  total 
vertical  stress  increment  can  be  given  by 


> 

(D 

1 

'  ^  -  M  K 

^-^vD  \ 

L  „  U  I  ""m 

<  ”p  ^p  / 

Substituting  Equations  6-6  and  6-7  into  Equation  6-15, 
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Let 


wL_  L_  _Ek^ 

ir 

g 


(6-17) 


fi  «  =  K  “ 

S  S  K  s 

g 


(6-18) 


Substituting  Equations  6-17  and  6-18  into  Equation  6-15  gives 


e  _  s  s 

^°VD  (fi  ^  +  K  )  +  K  ) 

s  ms  m 


Defining  the  following  nondimensional  quantities; 


(6-19) 


(6-20) 


(6-21) 


Equation  6-19  can  be  expressed  in  the  following  form 


1 

» 

< 

(1  - 

r) 

(1  +  a) 

(1  +  |) 

(6-22) 
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Partially  Coupled  Model  With  Constant  Poisson's  Ratio 

As  a  special  case,  assume  Poisson's  ratio  (y)  remains 
constant  during  the  dynamic  loading  and  unloading.  Then, 


Substitution  of 


K 


K 


s 


s 


(1  +  V) 
3(1  -  u) 


(1  u) 
3' a  -  U) 


Equations  6-23  and  6-24 


M  ^  (6-23) 

s 

M  ”  (6-24) 

s 

into  Equation  6-20  gives 

=  r  (6-25) 


That  is,  r  is  equal  to  the  strain  recovery  ratio,  r,  in  a  one 
phase  material.  Substitution  of  Equation  6-23  into  Equation  6-21 
gives 


(1  +  y) 

3(1  -  y)  Kg/ 


(6-26) 


Let 


a  = 


(6-27) 
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n  =  1 


(1  -»•  u) 

3(1  -  u)  Kg 


(6-28) 


Then 


a  =  a*n  (6-29) 


Substituting  Equations  6-26  and  6-29  into  Equation  6-22  gives 


AiTe 

(1  -  r) 

(1  +  an)  (1  + 

an 

Decoupled  Model 


(6-30) 


Residual  pore  pressure  response  can  be  obtained  in  a 
similar  manner  for  the  simpler  decoupled  model  described  by 
Blouin  and  Kim  (1984).  In  the  decoupled  model  the  compressibility 
of  the  solid  grains  by  the  pore  water  is  neglected.  The 
undrained  decoupled  loading  and  unloading  constrained  moduli  are 
given  by 


+  K 


m 


(6-31) 
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u 


(6-32) 


+  K 


m 


Using  the  same  procedure  as  employed  in  the  partially  coupled 
case,  the  following  stress  ratio  is  obtained  for  the  decoupled 
case: 


ATTe 

(1  - 

r) 

(1  +  a) 

(1  +  1  ) 

(6-33) 


Equation  6-33  has  the  same  form  as  Equation  6-22  except  that  r  and 
a  in  Equation  6-22  are  replaced  by  r  and  a  respectively. 


Post  Shot  Effective  Stress  and  Pore  Pressure 

Once  residual  excess  pore  pressure  (AiTg)  is  obtained,  the 
pore  pressure  ( rr)  and  effective  vertical  stress  (cj^')  at  the 
end  of  the  dynamic  unloading  are  calculated  as  follows: 

TT  =  ^  (6-34) 


a 


V 


(6-35) 


where  and  at’©  the  pre-shot  in  situ  pore  pressure 

and  effective  vertical  stress  respectively.  It  should  be  noted 
that  if  the  post-shot  effective  vertical  stress  (a^')  reaches 
the  tensile  strength  of  the  material,  the  material  would  be 
liquefied.  In  such  a  liquefied  material,  the  effective  stress  is 
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zero  and  the  post-shot  pore  pressure  (it)  reaches  a  maximum  which 
is  equal  to  the  in  situ  total  vertical  stress  (0^^). 

Parametric  Analysis 

A  parametric  study  was  conducted  to  examine  the  influence 
of  the  strain  recovery  ratio  (r)  and  the  ratio  of  the  skeleton 
loading  constrained  modulus  to  the  mixture  modulus  (a)  on  the 
residual  excess  pore  pressure  (Aitg). 

Equation  6-30,  which  represents  the  partially  decoupled 
model  with  constant  Poisson's  ratio,  was  used  in  this  parameter 
study.  Typical  material  parameters  representative  of  Generic 
Mountain  C  are  summarized  below. 

Porosity 
Poisson's  ratio 
Bulk  Modulus 

Pore  Water 

Solid  Grain 

Mixture 
n  =  0.8823 

Figure  6.5  summarizes  the  results  of  the  parametric 
study.  Residual  excess  pore  pressure  normalized  by  the  dynamic 
total  vertical  stress  increment  is  plotted  as  a  function  of 
strain  recovery  ratio  at  5  different  values  of  a  ranging  from  0.5 
to  2.5.  For  the  case  where  the  skeleton  loading  constrained 
modulus  (Mg^)  is  equal  to  the  mixture  modulus 
the  residual  excess  pore  pressure  is  20%  of  the  dynamic  total 
vertical  stress  increment  for  a  strain  recovery  ratio  of  0.46. 


n  =  20% 
y  =  0.2 

=  0.29  X  10®  psi 
Kg  =  5  X  10®  psi 
Kjji  =  1.177  X  10®  psi 
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The  normalized  excess  pore  pressure  is  10%  for  a  recovery  ratio 
of  0.67. 
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Permanent  Dispiacement  (ft) 


Figure  6.1.  Comparison  of  Rainier  radial  permanent 

displacements  to  theoretical  approximation. 
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Figure  6.2.  Residual  stress  fields  for  nearly  saturated  tuff 
events  normalized  by  the  final  cavity  radius  size 

(from  Patch,  1984). 
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Vertical  Stress  (6^ ) 


=  In  Situ  Effective  Vertical  Stress 

=  Loading  Constrained  Modulus  of 
Solid  Skeleton 

=  Unloading  Constrained  Modulus  of 
Solid  Skeleton 


Figure  6.3.  Bilinear  stress-strain  curve  for  solid  skeleton 
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Total  Vertical  Stress 


a  .  =  In  Situ  Total  Vertical  Stress 

VI 

Ao^D  *  Dynamic  Total  Vertical  Stress  Increment 
Ae^^  =  Volumetric  Strain  at  the  end  of  Loading 
Ae  ^  =  Recovered  Volumetric  Strain  at  the  end  of 

V 

Unloading 


-  Loading  Undrained  Constrained  Modulus  of  Bulk 
Mixture 


=  Unloading  Undrained  Constrained  Modulus  of 
Bulk  Mixture 


Figure  6.4.  Bilinear  stress-strain  curve  for  partially 
coupled  xmdrained  bulk  mixture. 
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Figure  6.5,  Normalized  residual  excess  pore  pressure  as  a  function  of  strain 
recovery  ratio. 


SECTION  7 


PARAMETRIC  EVALUATION  OF  FLOW  AND  PRESSURE  DISSIPATION 


CALCULATIONAL  AND  MATERIAL  PARAMETERS 


A  matrix  of  calculations  was  performed  to  determine 
potential  flow  into  tunnels  from  residual  loadings  similar  to 
those  developed  from  underground  nuclear  detonations  and  similar 
to  those  which  might  be  expected  from  surface  or  shallow  buried 
nuclear  detonations.  Material  properties  of  the  host  rock  were 
varied  over  a  broad  range  which  encompassed  nearly  the  entire 
range  of  properties  listed  for  Generic  Mountain  C  as  well  as  the 
range  of  properties  of  Rainier  Mesa  Tuff.  These  calculations  are 
primarily  meant  to  identify  potential  flow  problems  and  the 
magnitude  of  these  problems  in  a  broad  range  of  geologies;  they 
are  not  meant  to  analyze  the  response  of  a  specific  target  in  a 
specific  geology. 

In  each  of  the  parametric  calculations  a  segment  of 
tunnel  in  fully  saturated  rock  is  subjected  to  a  uniform  long 
term  loading.  This  loading  either  involves  imposition  of  a  total 
stress  on  the  grid  boundary  or  imposition  of  a  fixed  boundary 
displacement.  The  former  is  similar  to  the  residual  stress 
loading  which  would  develop  from  the  dynamic  loading  of  a 
hysteretic  saturated  rock  mass,  as  described  in  Section  5. 
Imposition  of  a  fixed  displacement  is  akin  to  the  loading  from 
an  underground  test,  where  the  rock  mass  is  subjected  to  a  fixed 
displacement  at  the  cavity  boundary.  Since  these  are 
axisymmetric  calculations,  however,  they  cannot  accurately 
replicate  the  loading  from  the  underground  test  conditions. 
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Figure  7.1  shows  a  schematic  section  view  of  the 
axisymmetr ic  loading  conditions  imposed  on  the  tunnel  throughout 
this  series  of  calculations.  A  20  ft  diameter  tunnel  was  used  in 
all  calculations.  Water  is  assumed  to  drain  freely  from  the 
tunnel  and  does  not  build  up  therein.  The  tunnel  walls  are  free 
to  displace  inward  under  the  imposed  loads  and  there  is  assumed 
to  be  no  lining  or  grout  in  the  vicinity  of  the  tunnel  to  inhibit 
flow  into  the  tunnel. 

In  all  but  the  initial  calculations,  the  loading  boundary 
is  assumed  to  be  impermeable  and  free  to  displace.  In  the 
initial  calculations,  a  fixed  displacement  was  imposed  on  the 
boundary  which  was  held  throughout  the  entire  time  span  of  the 
calculations.  As  shown  in  Figure  7.2,  flow  was  slightly  lower  in 
the  fixed  boundary  calculation.  This  is  expected  because  the 
boundary  does  not  move  inward  during  the  calculation  to 
accommodate  the  consolidation  near  the  tunnel  wall.  Thus,  pore 
pressures  are  somewhat  lower  and  total  stresses  dissipate  during 
the  fixed  boundary  calculation.  In  the  fixed  boundary 
calculation  the  imposed  displacement  is  equal  to  the  initial 
displacement  under  the  1000  psi  loading.  Thus,  effective 
stresses  and  pore  pressures  are  initially  the  same  in  both 
calculations.  Even  though  the  loading  boundary  is  free  to 
displace  inward  during  the  fixed  stress  calculation,  there  was  so 
little  difference  between  the  two  calculations  that  for 
convenience  the  fixed  displacement  loadings  were  not  used  in  this 
initial  study. 

The  range  to  the  loaded  boundary  was  varied  as  necessary 
in  each  calculation  to  insure  that  pore  pressure  at  the  boundary 
did  not  decrease  significantly  during  the  30  day  total  time  span 
of  the  calculation.  In  other  words,  in  more  permeable  material. 
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where  pore  pressures  dissipated  more  rapidly  and  further  into  the 
host  rock,  the  loaded  boundary  was  located  at  a  greater  range. 

The  range  to  the  loaded  boundary  was  adjusted  so  that  at  30  days 
the  pore  pressure  drop  at  the  boundary  was  less  than  5%  of 
the  initial  pore  pressure.  A  boundary  radius  of  7000  ft  was  used 
in  most  calculations,  though  this  was  increased  to  accommodate  the 
higher  permeability  materials. 

The  number  of  grid  elements  varied  in  proportion  to  the 
distance  to  the  loaded  boundary.  For  most  calculations,  there 
were  120  elements  between  the  tunnel  wall  and  the  loaded 
boundary.  However,  as  many  as  400  elements  were  used  in  the 
highest  permeability  calculations.  The  element  sizes  increased 
with  increasing  range  from  the  tunnel.  In  all  cases  half  of  the 
elements  were  within  one  quarter  of  the  distance  from  the  tunnel 
to  the  loaded  boundary.  The  element  closest  to  the  tunnel  was 
0.83/n  percent  of  the  total  grid  length  (where  n  is  the  total 
number  of  elements),  with  each  successive  element  increasing  in 
length.  Thus,  for  the  standard  120  element  7000  ft  grid,  the 
first  element  was  0.49  ft  long  and  there  were  60  elements  making 
up  the  initial  1750  ft  of  material  adjacent  to  the  tunnel. 

A  total  stress  of  1000  psi  was  applied  in  all 
calculations  of  the  parametric  matrix.  This  is  representative  of 
residual  stresses  in  the  vicinity  of  tunnels  at  ranges  of 
interest  in  the  underground  nuclear  shots.  It  is  also  a 
reasonable  approximation  of  the  residual  stresses  which  might 
develop  around  a  deep  based  system  from  nonrecover able 
deformation  of  the  surrounding  rock  due  to  surface  burst  loadings 
as  discussed  in  Section  5.  A  series  of  check  runs  was  performed 
in  which  the  applied  stress  was  varied  between  100  and  10,000 
psi.  For  all  material  properties  used  in  this  parametric 
analysis,  the  flow  and  pressure  dissipation  were  found  to  be  in 
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direct  proportion  to  the  applied  stress.  For  the  100  psi 
loading,  flow  and  pressures  were  one  tenth  those  in  the  1000  psi 
loading,  while  flows  and  pressures  were  ten  times  as  large  in  the 
10,000  psi  calculation.  Thus,  flows  and  pressures  for  loading 
stresses  other  than  1000  psi  can  be  easily  calculated  by 
multiplying  the  calculated  results  bv  ratio  between  the 

desired  applied  stress  and  the  lOOu  Dad. 

A  suite  of  rock  properties  encompassing  nearly  the  entire 
range  of  properties  found  in  the  Rainier  Me-a  tuffs  and  in  the 
various  rocks  of  Generic  Mountain  C  were  incluled  'n  tl»e 
parametric  calculations.  For  purposes  of  this  initial  stt;Qy, 
both  the  rock  skeleton  and  pore  water  were  assumed  to  be  linear 
and  elastic.  Material  properties  used  in  the  calculations  are 
summarized  in  Table  7.1.  Properties  held  constant  throughout  all 
the  various  parameter  studies  included  the  100%  degree  of 
saturation,  the  bulk  modulus  of  water  with  a  value  of  0.29  x 
10®  psi,  and  the  bulk  modulus  of  the  solid  grains  with  a 
value  of  5  X  10®  psi.  Material  parameters  which  were  varied 
in  the  calculations  include  the  permeability,  Poisson's  ratio, 
constrained  modulus  and  porosity  of  the  rock  skeleton. 

Values  of  permeability  had  the  widest  range  of  any 
parameter  investigated  in  this  study.  Permeability  ranged  from  a 
low  value  of  6.7  x  10”®  ft/day  (9.3  x  10“^®  in/s), 
representative  of  the  lower  values  reported  for  the  zeolitized 
Rainier  Mesa  tuffs,  to  100  ft/day  (1.39  x  10”^  in/s) 
representative  of  the  very  permeable  Generic  Mountain  C 
sandstones.  Intermediate  values  represent  more  permeable  or 
fractured  tuffs  and  the  various  other  formations  of  Mountain  C. 

A  standard  intermediate  permeability  of  0.1  ft/day  (1.39  x 
10”®  in/s)  was  used  in  the  parametric  evaluations  of  the 
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influence  of  Poisson's  ratio/  porosity,  and  skeleton  modulus  on 
flow  and  pore  pressure  dissipation. 

The  constrained  modulus  of  the  rock  skeleton  was  varied 
between  0,5  x  10®  and  5  x  10®  psi.  The  lowest  value  is 
representative  of  soft  altered  tuff  while  the  highest  is 
representative  of  the  hardest  most  competent  rocks  of  Mountain  C. 

A  standard  value  of  2.5  x  10®  psi  was  used  in  the  studies  of 
variation  in  porosity  and  Poisson's  ratio,  while  a  value  of  1.9  x 
10®  psi  was  used  in  the  permeability  study.  Either  value  is  a 
reasonable  representation  of  the  average  properties  of  Generic 
Mountain  C. 

Poisson's  ratio  of  the  rock  skeleton  was  varied  from  0.1 
to  0,4,  with  a  standard  value  of  0.2  used  in  the  other  parameter 
studies.  The  lower  value  represents  the  harder  more  competent 
Mountain  C  rocks  while  the  upper  limit  approximates  that  in  the 
weaker  zeolitized  tuffs.  The  standard  value  is  representative  of 
Generic  Mountain  C. 

The  final  parameter  varied  in  this  study  was  the  skeleton 
porosity.  A  value  of  0.2  was  used  as  the  standard.  This  is 
representative  of  average  values  for  Generic  Mountain  C.  Other 
values  ranged  from  0.05,  representative  of  dense  well  consolidated 
sedimentary  rock  to  0.4,  representative  of  very  porous  tuff. 

The  influences  of  the  above  parameter  variations  on  flow 
and  pore  pressure  dissipation  are  discussed  in  the  following 
subsections . 
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INFLUENCE  OF  PERMEABILITY 


The  influence  of  permeability  on  cumulative  flow  and  flow 
rate  under  the  1000  psi  total  stress  loading  is  summarized  in 
Figures  7.3  and  7,4.  These  figures  show  cumulative  flow  and  flow 
rate  per  foot  of  tunnel  as  a  function  of  permeability  at  fixed 
times  of  1  hr  and  4  hr,  and  1,  7  and  30  days.  The  most  striking 
feature  of  these  plots  is  the  large  flows  and  flow  rates  which 
develop  in  the  medium  to  high  permeability  rock.  Topi  et  al. 
(1984)  indicate  a  typical  range  of  permeabilities  for  rock  in 
Generic  Mountain  C  as  0.5  to  50  ft/day  (6.944  x  10“^  to  6.944 
X  10“^  in/s).  Initial  flows  during  the  first  day  of  between 
5000  gallons  and  400,000  gallons  per  foot  of  tunnel  would  be 
expected  for  this  permeability  range.  Since  it  takes  only  2350 
gallons  per  foot  to  completely  fill  the  tunnel,  under  these 
assumed  conditions  there  appears  to  be  a  very  serious  potential 
flow  problem  for  a  deep  based  system  at  such  a  site.  This  range 
in  cumulative  flows  equates  to  an  average  first  day  flow  rate  of 
from  3.5  to  278  gpm  per  foot  of  tunnel. 

Contrast  these  values  to  the  one  day  flows  expected  in 
the  saturated  zeolitized  tuffs.  In  the  tuff,  first  day  flows  of 
between  1.2  to  10  gallons  per  foot  of  tunnel  would  be  expected 
for  the  typical  range  of  tuff  permeabilities  of  Table  2.1.  These 
equate  to  flow  rates  of  only  8  x  10”^  to  7  x  10“^  gpm  per 
foot  of  tunnel.  At  the  lower  end  of  the  tuff  permeabilities,  the 
incoming  water  would  likely  evaporate  as  fast  as  it  entered  the 
tunnel,  and  toward  the  upper  end  of  the  permeability  range  the 
water  would  probably  pose  only  a  slight  inconvenience  rather  than 
a  problem.  Minimal  inflow,  as  indicated  by  these  calculations, 
has  been  experienced  in  Rainier  Mesa,  except  where  fracture  water 
has  been  mobilized.  The  only  measurements  of  nonfracture  inflow. 
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i.e.  squeezing  of  interstitial  water  from  the  rock  pores,  were 
made  on  ONETON.  These  measurements  compare  favorably  with 
calculated  flows  as  demonstrated  in  Section  8. 

Examples  of  calculated  flow,  flow  rate  and  pore  pressure 
dissipation  in  the  lowest  permeability  rock  examined  in  this 
parameter  study  (k  =  9,3  x  10“^^  in/s)  are  shown  in  Figures 
7.5  through  7.7,  Accumulated  flow  and  flow  rate  per  foot  of 
tunnel  are  shown  as  functions  of  time  for  the  full  30  days  in  the 
first  two  figures,  while  pore  pressures  as  a  function  of  range 
are  shown  at  various  times  from  10  sec  to  30  days  in  Figure  7.7. 

In  this  low  permeability  rock  only  about  12  gallons  of  water 
accumulate  per  foot  of  tunnel  over  the  entire  30  day  time  span. 

The  relatively  high  initial  flow  rates  decrease  very  rapidly 
within  the  first  several  hours  to  a  very  low,  slowly  decaying  rate 
during  the  bulk  of  the  time  span.  Pore  pressure  dissipation 
advances  very  slowly  into  the  host  rock.  At  30  days  pore 
pressures  have  begun  to  decay  only  to  a  range  of  65  ft  from  the 
tunnel  wall. 

Similar  plots  for  an  intermediate  permeability  rock  (k  = 
1.39  X  10”^  in/s)  are  shown  in  Figures  7.8  through  7.10. 

This  intermediate  permeability  is  near  the  lower  bound  value 
assumed  for  Generic  Mountain  C.  Even  in  this  intermediate 
permeability  rock,  flows  and  flow  rates  are  reaching  rather 
alarming  proportions  in  terms  of  a  deep  based  facility. 

Accumulated  flow  exceeds  30,000  gallons  per  foot  of  tunnel  over 
the  30  day  time  span.  Flow  rate  is  initially  about  3  gpm  per 
foot  of  tunnel,  but  drops  only  modestly  to  about  0.6  gpm  at  the 
end  of  30  days.  The  pore  pressure  dissipation  front  advances 
very  rapidly  into  the  host  rock,  reaching  a  range  of  about  1600 
ft  from  the  tunnel  at  the  end  of  the  first  day  and  6000  ft  at  the 
end  of  30  days.  This  rapid  dissipation  suggests  that  finite 
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geologic  and  geometric  parameters  in  an  actual  deep  based  system 
will  themselves  limit  flow  into  the  tunnels;  i.e.  factors  such  as 
finite  depths  of  cover  and  limited  aquifers  will  somewhat 
restrict  the  high  accumulated  flows  and  flow  rates  indicated  at 
later  times  in  these  calculations.  Boundary  limitations  which 
might  replicate  such  restrictions  were  not  used  in  this  analysis 
so  as  not  to  confuse  the  comparisons  with  the  more  impermeable 
rocks. 


Figures  7.11  through  7.13  present  the  accumulated  flows, 
flow  rates  and  pore  pressure  dissipations  in  the  highest 
permeability  rock  (k  =  1.39  x  10  in/s).  Exceptionally  high 
flows  and  flow  rates  occur  in  this  highly  permeable  rock,  with 
pore  pressure  dissipation  extending  for  tens  of  thousands  of 
feet.  Again,  these  calculations  do  not  have  realistic  boundary 
conditions,  so  the  large  flows  indicate  only  that  a  substantial 
problem  exists.  They  should  not  be  used  to  estimate  actual 
expected  flow  magnitudes. 


INFLUENCE  OF  SKELETON  MODULUS 

The  influence  of  varying  skeleton  constrained  modulus 
between  0.5  x  10®  psi  and  5  x  10®  psi  on  accumulated 
flow,  flow  rate  and  pore  pressure  dissipation  is  shown  in  Figures 
7.14  through  7.16.  In  this  set  of  calculations,  the  intermediate 
permeability  of  1.39  x  10“®  in/s  was  used.  Other  material 
properties  are  listed  in  Table  7.1.  As  would  be  expected,  in 
the  stiffer  rocks,  the  skeleton  assumes  a  greater  percentage  of 
the  imposed  load.  Thus,  the  generated  pore  pressures  and 
resultant  flows  are  lower  in  the  stiffer  rocks.  Accumulated  flow 
at  30  days  in  the  softest  rock  is  about  4.7  times  greater  than 
flow  in  the  stiffest  rock.  The  large  differences^  in  initial  pore 
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pressure  are  evident  in  Figure  7.i6a,  where  the  initial  pore 
pressure  in  the  softest  rock  is  about  4.3  times  that  in  the 
hardest  rock.  These  initial  pore  pressures  can  be  calculated 
directly  from  Equation  A-54  for  the  undrained  loading  condition. 
Figure  7.17  is  a  plot  of  pore  pressure  as  a  function  of  skeleton 
modulus  for  an  undrained  loading  of  1000  psi.  This  plot  was 
generated  from  Equation  A-54  using  the  material  properties  from 
Table  7.1. 

INFLUENCE  OF  POISSON'S  RATIO 

Poisson's  ratio  of  the  rock  skeleton  was  found  to  have 
only  a  modest  influence  on  flow  and  flow  rate.  Figures  7.18 
through  7.20  show  accumulated  flow,  flow  rate  and  pore  pressure 
dissipation  for  skeletons  having  Poisson's  ratios  ranging  from 
0.1  to  0.4.  As  indicated  in  Table  7.1,  the  constrained  moduli 
and  permeabilities  were  assumed  constant  at  2.5  x  10'^  and 
1.39  X  10“^  in/s  respectively  in  these  calculations.  The 
accumulated  flow  at  30  days  in  the  rock  with  a  Poisson's  ratio  of 
0.1  is  about  50%  greater  than  flow  in  the  rock  having  a  Poisson's 
ratio  of  0.4.  The  corresponding  flow  rates  and  initial  pore 
pressures  are  also  higher  in  the  low  Poisson's  ratio  rock.  The 
initial  undrained  pore  pressures  can  be  calculated  directly  from 
Equation  A-54.  Note,  however,  that  Equation  A-54  uses  the 
skeleton  bulk  modulus  which  is  a  function  of  the  constrained 
modulus  and  Poisson's  ratio  according  to 


K  =  M  (1  +  y) 

®  ®  3(1  -  y) 
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Thus/  for  a  constant  constrained  modulus,  the  bulk  modulus 
increases  with  increasing  Poisson's  ratio.  The  combined  effect 
of  the  increasing  bulk  modulus  and  increasing  Poisson's  ratio  is 
a  moderate  decrease  in  initial  pore  pressure  with  increasing 
Poisson's  ratio.  This  effect  is  shown  in  Figure  7.21  where  pore 
pressure  is  plotted  as  a  function  of  Poisson's  ratio  for  a 
constant  value  of  constrained  modulus.  In  contrast,  for  a 
constant  skeleton  bulk  modulus,  initial  pore  pressure  increases 
slightly  with  increasing  Poisson's  ratio.  The  initial  pore 
pressure  as  a  function  of  Poisson's  ratio  for  a  constant  skeleton 
bulk  modulus  of  1.25  x  10®  psi  is  plotted  for  comparison  in 
Figure  7.21. 

INFLUENCE  OF  POROSITY 


The  final  property  examined  in  this  parameter  study  was 
the  influence  of  porosity  on  flow.  These  results  gave 
relationships  which  seemed  opposite  to  those  one  would 
intuitively  expect.  As  shown  in  Table  7.1,  porosity  was  varied 
from  5%  to  50%  while  permeability,  skeleton  constrained  modulus 
and  Poisson's  ratio  were  held  constant  at  values  of  1.39  x 
10“^  in/s,  2.5  X  10®  psi  and  0.2  respectively.  The 
resulting  plots  of  accumulated  flow,  flow  rate  and  pore  pressure 
dissipation  in  Figures  7.22  through  7.24  show  that  pore  pressures 
and  the  resulting  flows  are  highest  when  porosities  are  lowest. 
This  nonintuitive  behavior  arises  because  porosity  is  varied 
independently  of  the  other  material  parameters.  Realistically, 
this  could  not  occur  over  such  a  great  a  range  in  porosity. 

The  initial  pore  pressure  response  as  a  function  of 
porosity  is  governed  by  the  Wood  equation  (Equation  A-13).  As 
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pcrosity  increases,  the  bulk  modulus  of  the  solid-water  mixture 
decreases  because  there  are  less  very  stiff  solid  grains 
in  the  mixture.  The  pore  pressure  response  of  Equation  A-54  is 
in  turn  directly  proportional  to  the  mixture  modulus,  i.e.  the 
higher  the  porosity,  the  lower  the  mixture  modulus  and  the  lower 
the  resultant  pore  pressure  under  the  undrained  loading. 
Therefore,  with  all  other  properties  held  constant,  the  pore 
pressures  and  flows  are  highest  when  permeabilities  are  lowest. 
The  accumulated  flow  at  30  days  in  the  lowest  porosity  rock  is 
2.9  times  that  in  the  highest  porosity  rock.  As  shown  in  figure 
7.24a,  this  is  because  the  initial  pore  pressure  in  the  low 
porosity  rock  is  about  3  times  that  in  the  high  porosity  rock. 
The  influence  of  porosity  on  the  initial  pore  pressure  is  shown 
in  Figure  7.25.  This  plot  was  computed  from  Equation  A-54  with 
material  constants  taken  from  Table  7.1. 


SUMMARY 


From  the  series  of  axisymmetric  calculations  used  to 
study  the  influence  of  material  property  variations  on  flow  and 
pore  pressure  dissipation  in  this  section  we  conclude  the  following; 

•  Flow  of  interstitial  pore  water  generated  by  explosively 
induced  residual  stresses  in  permeable  saturated  or  nearly 
saturated  rocks  such  as  those  in  Generic  Mountain  C  is  a 
serious  potential  threat  to  deep  based  systems; 

•  Interstitial  flow  generated  by  similar  loadings  in  non- 
permeable  rock  such  as  the  zeolitized  tuffs  of  Rainier 
Mesa  is  minimal  and  both  experience  and  calculations 
indicate  that  this  type  of  flow  would  not  be  a  major 
problem  to  a  deep  based  system  in  similar  rock; 
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•  Interstitial  pore  water  flow  is  extremely  sensitive 
to  large  changes  in  permeability/  is  somewhat  sensi¬ 
tive  to  skeleton  modulus  and  porosity,  and  is  rel¬ 
atively  insensitive  to  skeleton  Poisson's  ratio  and 
to  whether  residual  stresses  are  constant  or  result 
from  a  fixed  boundary  displacement. 
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Table  7.1.  Material  properties  used  in  parametric  calculations. 


PROPERTY 

STANDARD  VALUE 

RANGE  OF  VALUES 

Skeleton  Constrained  Modulus 

2.56  X  10^  psi  or 

0-56  X  10^  to 

1.9  X  10^  psi 

5.0  X  10^  psi 

Bulk  Modulus  of  Pore  Water 

0.29  X  10^  psi 

Skeleton  Porosity 

0.2 

0.05  to  0.5 

Skeleton  Permeability 

0.1  ft/day 

6.7  X  10  ^  to 
100  ft/day 

1.39  X  10  ^  in/s 

9.28  X  10~^°  to 
1.39  X  10  ^  in/i 

Degree  of  Saturation 

100% 

Skeleton  Poisson's  Ratio 

0.2 

0.1  to  O.A 

Bulk  Modulus  of  Solid  Grains 

5  X  10^  psi 
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Figure  7.2.  Comparison  of  cumulative  stress  from  total  stress  and  fixed  displacement 
loadings. 
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Figure  7.4.  Flow  rate  as  a  function  of  permeability. 
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.285E-10  in/s 


Figure  7.7a.  Pore  pressure  dissipation  in  low  permeability  rock 
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Cumulative  flow  in  inedium  permeability  rock 
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Figure  7.10a.  Pore  pressure  dissipation  in  medium  permeability  rock 
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Figure  7.11.  Cumulative  flow  in  high  permeability  rock. 
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Figure  7.13a.  Pore  pressure  dissipation  in  high  permeability  rock 


>000.  50000.  75000.  100000.  125000.  150000 


Figure  7.16a.  Influence  of  skeleton  modulus  on  pore  pressure 
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Range  (ft) 


Pore  Pressure 


Poisson's  Ratio 

Figure  7.21.  Initial  pore  pressure  as  a  function  of 
Poisson's  ratio. 
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Figure  7.22.  Influence  of  porosity  on  flow  volume 
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Figure  7.25.  Initial  pore  pressure  as  a  function  of 
porosity. 


SECTION  8 


ONETON  NUMERICAL  SIMULATION 


INTRODUCTION 

Calculations  in  Section  7  indicated  that  flow  of 
interstitial  pore  water  due  to  blast  induced  residual  stresses  in 
the  permeable  rocks  of  Generic  Mountain  C  poses  a  serious 
potential  threat  to  deep  based  systems.  However,  in  the  much 
less  permeable  tuffs  of  Rainier  mesa,  flow  of  interstitial  pore 
water  has  only  been  a  minor  problem.  Pore  water  flow  in  the 
Rainier  Mesa  tuff  has  only  been  well  documented  on  one  event,  the 
ONETON  HE  experiment,  the  data  from  which  is  summarized  in 
Section  4.  In  order  to  give  credibility  to  the  Generic  Mountain 
C  calculations,  the  CONSL  axisymmetric  code  was  used  to  calculate 
the  flow  of  water  on  the  ONETON  event.  If  the  ONETON  results 
could  be  replicated  using  measured  residual  stresses  and  actual 
tuff  properties,  it  was  felt  that  the  trends  toward  very  high 
flows  in  the  Generic  Mountain  C  calculations  become  much  more 
credible.  The  results  of  the  ONETON  calculations  are  presented 
in  this  section. 

PROBLEM  CHARACTERIZATION  AND  RESULTS 

Since  the  CONSL  code  is  only  a  one  dimensional 
axisymmetric  code,  some  simplifying  assumptions  with  regard  to 
the  problem  geometry  and  residual  stress  loading  had  to  be  made. 
Plan  and  section  views  of  the  idealized  problem  geometry  are 
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shown  in  Figure  8.1.  The  drift  into  which  the  ONETON 
instrumentation  cables  were  strung,  pictured  in  Figure  4.6,  was 
approximated  as  a  vertical  right  circular  cylinder.  The  16  ft 
diameter  and  8  ft  height  approximated  the  actual  dimensions  of  the 
drift.  The  center  of  the  drift  was  at  a  range  of  63  ft  from  the 
working  point,  which  placed  the  closest  face  at  the  actual  range 
of  55  ft. 

A  symmetric  total  residual  stress  of  630  psi  was  applied 
to  the  remote  grid  boundary.  This  approximated  the  residual 
stress  at  the  63  ft  range  of  the  drift  centroid.  The  630  psi 
residual  stress  loading  was  obtained  from  a  linear  fit  to,  and 
extrapolation  of,  the  HE  residual  stress  data  presented  in  Figure 
4.8.  This  data  fit  is  shown  in  Figure  8.2.  In  the  actual  case 
the  loading  would  be  considerably  more  complex,  with  stresses  on 
the  working  point  side  being  considerably  in  excess  of  those 
assumed  here  and  with  stresses  on  the  far  side  of  the  drift  being 
considerably  less. 

The  ONETON  calculations!  parameters  and  material 
properties  used  in  the  CONSL  calculation  are  listed  in  Table  8.1. 

A  calculational  grid  containing  60  elements  extended  to  a  range 
of  200  ft  from  the  tunnel  wall.  The  630  psi  total  stress  was 
applied  to  this  impermeable  remote  boundary.  The  material 
properties  of  the  tuff  were  typical  zeolitized  tuff  properties 
extracted  from  a  report  by  Blouin  and  Kim  (1983).  The 
permeability  was  iteratively  varied  until  the  pore  water  flow 
matched  the  measured  ONETON  flows  reported  by  Smith  (1983).  Flow 
per  unit  area  was  computed  for  a  circular  section  of  drift  and 
was  then  multiplied  by  the  total  area  of  the  idealized 
cylindrical  drift  to  give  the  total  accumulated  flow.  The 
permeability  of  2.8  x  10”®  in/s,  which  gave  the  best  match  to 
the  ONETON  flow  data,  was  near  the  upper  end  of  the  range  of  the 
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tuff  permeabilities  reported  by  Thordarson  (1965)  and  shown  in 
Figure  7.3. 


The  calculated  cumulative  flow  and  pore  pressure 
dissipation  profiles  for  the  ONETON  event  are  shown  in  Figures 
8.3  and  8.4  respectively.  The  flows  reported  by  Smith  are 
superimposed  on  the  calculation.  The  first  flow  measurement  is 
somewhat  low  because  some  drying  out  of  the  pore  water  near  the 
drift  walls  probably  occurred  prior  to  the  detonation.  This  was 
not  modeled  in  the  calculation.  No  additional  flow  measurements 
were  obtained  followi  ig  the  second  measurement  10  days  after  the 
event. 


Even  though  rather  crude  simplifying  assumptions  had  to 
be  made  in  order  to  calculate  the  ONETON  flow  using  the  CONSL  one 
dimensional  code,  the  calculated  flows  were  a  good  match  to  the 
actual  flows.  The  fact  that  this  agreement  was  achieved  using 
material  properties  that  are  within  the  best  estimates  of  the 
tuff  properties  and  applied  loads  which  were  obtained  from  test 
data  is  strong  validation  of  the  calculat ional  procedure.  We 
believe  this  agreement  supports  Smith's  hypothesis  that  ONETON 
pore  water  flow  was  caused  by  the  residual  stress  field  forcing 
the  pore  water  from  the  voids.  The  agreement  between  the 
calculations  and  the  field  data  gives  real  credibility  to  the 
conclusion  that  residual  stress  induced  flow  is  a  serious 
potential  problem  for  deep  based  sites  in  permeable  rocks  such 
those  of  Generic  Mountain  C. 
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Table  8.1.  ONETON  calculational  parameters  and  material 
properties . 


Tunnel  Radius 

Radius  of  Remote  Boundary 

Number  of  Elements 

Mesh  Growth  Factor 

Applied  Total  Stress 

Permeable  Tunnel  Boundary 

Impermeable  Remote  Boundary 


Skeleton  Constrained  Modulus 
Skeleton  Poisson's  Ratio 
Porosity 
Permeability 

Bulk  Modulus  of  Pore  Water 
Bulk  Modulus  of  Solid  Grains 


8.  ft 
208.  ft 


0.25 
630  psi 


1.07  X  10°  psi 

0.36 

0.34 

2.8  X  10~®  in/s 
0.29  X  10®  psi 
5.  X  10®  psi 
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Figure  8.1.  Schematic  views  of  idealized  ONETON  layout. 
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Figure  8.2.  Fit  to  HE  residual  stress  data. 
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Figure  8.4.  Calculated  CNETON  pore  pressure  dissipation  profiles 
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APPENDIX  A 


UNDRAINED  ISOTROPIC  COMPRESSION  OF  SATURATED  POROUS 

MEDIA  IN  PLANE  STRAIN 


INTRODUCTION 

Relationships  between  effective  stresses,  total  stresses 
and  pore  pressure  during  undrained  isotropic  and  constrained 
(uniaxial  strain)  loadings  of  saturated  porous  elastic  materials 
were  derived  by  Blouin  and  Kim  (1984).  These  relationships  are 
fundamental  to  the  development  of  sophisticated  multiphase 
computer  codes  used  to  model  the  behavior  of  saturated  soil  and 
rock.  Three  coupling  relationships  between  the  material  skeleton 
and  pore  water  were  postulated.  The  first,  termed  the  decoupled 
model,  is  the  simplest  and  generally  the  least  accurate  of  the 
three.  In  the  decoupled  model  the  compressibility  of  the 
material  skeleton  and  that  of  the  pore  water  are  assumed  to  act 
in  parallel,  but  completely  independent  of  one  another.  That  is, 
the  material  skeleton  is  compressed  solely  by  the  intergranular 
or  effective  stresses  and  the  mixture  of  pore  water  and  solid 
grains  is  compressed  solely  by  the  pore  pressure.  The  volume 
compression  of  each  phase  is  equal  to,  but  independent  of,  the 
compression  of  the  other  phase. 

In  the  second  more  sophisticated  model,  termed  the 
partially  coupled  model,  the  compressibility  of  the  skeleton  is 
linked  to  that  of  the  pore  water  in  that  the  pore  water  pressure 
compresses  the  solid  grains  comprising  the  skeleton.  This 
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results  in  a  contraction  of  the  skeleton  in  addition  to  the 
skeleton  volume  reduction  caused  by  the  effective  stresses. 

The  third  and  most  sophisticated  model,  termed  the  fully 
coupled  model,  utilizes  the  features  of  the  partially  coupled 
model,  but  also  further  links  the  skeleton  and  pore  water 
response  by  taking  into  account  the  compression  of  the  solid 
grains  by  the  effective  stresses  in  the  soil  skeleton.  This 
compression  is  in  addition  to  the  volume  change  due  to  the 
compression  of  the  solid  grains  and  pore  water  by  the  pore 
pressure.  The  latter  is  a  feature  of  the  solid  grain/pore  water 
mixture  compressibility  of  all  three  models. 

The  isotropic  and  uniaxial  loading  conditions  treated  by 
Blouin  and  Kim  (1984)  are  not  applicable  to  many  problems 
including  plane  strain  loading  conditions.  Plane  strain 
Conditions  are  often  used  in  the  analysis  of  long  structures  and 
loadings  such  as  those  around  tunnels  and  footings  where  there 
is  assumed  to  be  no  strain  in  the  axial  direction.  Also,  two 
dimensional  computer  codes  often  approximate  three  dimensional 
problems  using  plane  strain  loading  conditions.  Because  of  the 
plane  strain  restriction,  the  results  from  these  calculations  may 
differ  significantly  from  those  expected  in  the  actual  problems 
of  interest.  In  this  Appendix,  relationships  between  stresses, 
pore  pressures  and  volume  strains  are  derived  for  undrained 
isotropic  loading  conditions  in  plane  strain.  Sets  of  equations 
for  each  of  the  three  coupling  models  are  presented  along  with 
comparisons  to  the  corresponding  equations  for  pure  isotropic 
loadings.  These  equations  were  used  to  check  pore  pressures  and 
effective  stresses  under  the  initial  loading  conditions  in  all 
the  CONSL  calculations. 
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DECOUPLED  MODEL 


The  assumed  plane  strain  loading  conditions  and  notation 
are  shown  in  Figure  A.l.  An  element  of  saturated  porous  material 
is  loaded  in  plane  strain  by  a  symmetric  stress  denoted  by  Oj.. 
Since  strain  in  the  axial  direction  is  restricted  to  zero,  a 
resultant  total  stress  in  the  axial  direction,  o^,  develops  which 
satisfies  the  zero  strain  restriction.  Each  of  the  total 
stresses  is  composed  of  an  intergranular  or  effective  stress 
component,  denoted  by  and  Og',  and  a  pore  pressure  u. 

The  pore  pressure  is  hydrostatic  and  acts  in  all  directions 
within  the  element.  According  to  the  effective  stress  law  the 
total  stresses  equal  the  effective  stresses  plus  the  pore 
pressure  according  to 


a 


r 


0^’  +  u 


(A-1) 


o 


a 


o'  +  u 

a 


(A-2) 


The  total  volume  change  in  the  element  due  to  application 
of  the  total  stresses  can  be  expressed  individually  in  terms  of 
both  the  pore  pressure  u  and  the  effective  stresses,  Oj-'  and 
Ojj'.  An  increase  in  pore  pressure  compresses  both  the  water 
within  the  granular  matrix  and  the  grains  themselves.  The  strain 
in  the  water  due  to  application  of  the  pore  pressure  u  is  given  by 
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AV 
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V 


w 


u 


K 


w 


(A-3) 


where  is  the  total  volume  of  pore  water  in  the  element  and 
Kw  is  the  bulk  modulus  of  the  pore  water.  The  total  volume 
of  water  in  an  element  of  material  having  a  porosity  n  and  a 
volume  V  is 


V  =  nV  (A-4) 

w 

The  volume  change  of  the  water  within  the  element  is  obtained  by 
combining  A-3  and  A-4 


AV^  =  n 


u 


V 


(A-5) 


Assuming  a  unit  initial  volume,  volume  strain  and  volume  change 
are  both  expressed  as 


AV  =  n  (A-6) 

w  K 

w 

The  pore  pressure  subjects  each  of  the  grains  making  up 
the  material  skeleton  to  a  hydrostatic  pressure,  u.  Each  grain 
undergoes  a  volume  strain  of 


u 


A4 


(A-7) 


where  AVgg  and  are  the  volume  change  and  total  volume 

of  a  single  grain  and  Kg  is  the  bulk  modulus  of  the  solid 
grains.  The  total  volume  change  of  the  grains  within  the  element 
is  given  by 


AV^  =  (1  -  n)  V 


(A-8) 


The  quantity  (1  -  n)V  is  the  initial  volume  of  the  solid  grains. 
For  a  unit  initial  element  volume,  the  volume  change  and  volume 
strain  in  the  solid  grains  are  given  by 


(1  -  n)  ^ 

g 


(A-9) 


The  total  volume  change  in  the  saturated  element,  AV, 
must  equal  the  sum  of  the  volume  changes  in  the  pore  water  and 
the  solid  grains  according  to 


AV  =  AV  +  AV 
w  g 


(A-10) 


The  total  volume  change  in  an  element  of  unit  volume  can  also  be 
expressed  as 


AV 


(A-11) 


where  Kj^^  is  the  bulk  modulus  of  the  solid  grain/pore  water 
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mixture.  Substitution  of  Equations  A-6,  A-9  and  A-11  into 
Equation  A-lO  gives 


(1  -  n) 


u 


(A-:2) 


Solution  of  A-12  gives  the  bulk  modulus  of  the  solid  grain/pore 
water  mixture  as 


K  K 
g  w 


K  +  n(K  -  K  ,) 
w  gw 


(A-13) 


which  is  equivalent  to  the  derivation  by  Wood  (1930)  for  a  soil 
water  suspension. 

The  volume  change  can  also  be  expressed  in  terras  of  the 
effective  stresses  acting  on  the  material  skeleton  as  ^depicted  in 
Figure  A-2.  The  plane  strain  loading  conditions  are  simulated  by 
summing  a  true  isotropic  loading  under  the  effective  stress 
Cj- '  and  a  uniaxial  tensile  stress  loading  of  -  Oj. ' 

in  the  axial  direction.  The  volume  change  from  the  isotropic 
loading  on  an  element  of  unit  volume  is  given  by 


(A-14) 


where  Kg  is  the  bulk  modulus  of  the  skeleton.  Kg  is  the 
modulus  that  would  be  obtained  from  a  fully  drained  hydrostatic 
loading  of  the  material  element.  The  volume  change  of  the 
skeleton  under  the  axial  loading  is  given  by 


A6 


2u 


(A-15) 


where  Eg  is  the  Young's  modulus  of  the  skeleton  and  y  is 
Poisson's  ratio  of  the  skeleton.  The  first  terra  in  A-15  is  the 
axial  extension  of  the  element  and  the  second  term  represents  the 
sum  of  the  two  components  of  radial  and  circumferential 
contraction.  Equation  A-15  is  simplified  as 


a,  “  o  ’ 

AVa  =  — —  (1  -  2y)  (A-16) 

Using  the  elastic  relationship  between  bulk  modulus  and  Youna's 
modulus  of 


E  =  3K  (1  -  2y)  (A-17) 

s  s 

Equation  A-17  can  be  expressed  in  terms  of  bulk  modulus  as 


(A-18) 


The  total  volume  change  of  the  skeleton  is  the  sum  of  the  volume 
changes  due  to  the  isotropic  and  axial  loadings  given  by 
Equations  A-14  and  A-18. 
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if 


AV  »  AV.  +  AV 

X  CL 


or 


AV 


2o  ’  +  a  ' 
r  a 


IF 


s 


The  mean  effective  stress  Ojjj* ,  is  given  by 


2a  *  +  o  ' 
•  -  r  a 


m 


(A-19) 


(A-20) 


(A-21) 


and  Equation  A-20  in  terms  of  the  mean  effective  stress  is  simply 


AV  = 


(A-22) 


Since  the  axial  strain  is  zero,  the  volume  change  is  the  sum  of 
the  two  strain  components  under  the  isotropic  stress  components, 


AV  =  2e^ 


{A-23) 


where  is  the  resultant  strain  from  each  of  the  components 
of  effective  stress,  a^' , 

The  final  compatibility  equation  is  obtained  from  the 
zero  strain  requirement  in  the  axial  direction.  From  Figure  A. 2, 
the  axial  skeleton  strain,  Sg,  equals  one  third  of  the  volume 
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strain  under  the  isotropic  loading  plus  the  axial  strain  from  the 
axial  loading; 


c 


+ 


(A-24) 


With  axial  strain  set  to  zero  and  Eg  expressed  in  terms  of 
Kg  according  to  Equation  A-17,  the  axial  effective  stress  is 
obtained  as  a  function  of  the  radial  effective  stress  according 
to 


o 


f 

a 


(A-25) 


The  five  equations  representing  the  decoupled 
compatibility  conditions  for  the  isotropic  loading  in  plane 
strain  are  summarized  below. 


o 

r 


+  u 


(A-l) 


o 


a 


a  '  +  u 
a 


(A-2) 


AV 


(A-11) 


AV  = 


+ 


(A-20) 
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{A-25) 


o 


f 

a 


=  2yOj.‘ 


These  are  solved  for  the  unknown  effective  stresses  Oj-'  and 
Og ' ,  the  pore  pressure,  u,  volume  change  AV,  and  the  total 
axial  stress,  all  as  a  function  of  the  prescribed 

isotropic  total  stress, 

Equating  equations  A-11  and  A-20  gives 


u 


K 


2a  '  +  a  • 
r  a 


m 


3ir 


(A-26) 


Substitution  of  Equation  A-25  into  A-26  yields 

2K 

U  =  ^  Oj.'  (1  +  U)  (A-27) 

and  substitution  of  A-27  into  A-1  gives  the  isotropic  effective 
stress  as 


3K 

_  •  _  „  _ s _ 

2K  (1  +  u)  +  3K^ 
m  s 


(A-28) 


Further  substitution  of  Equation  A-28  into  Equation  A-1  gives  the 
pore  pressure  as 


AlO 


u 


(A-29) 


2Kjj^(l  +  p) 

°r  2K  (1  +  p)  +  3K 
in  s 


The  effective  axial  stress  is  obtained  by  substitution  of 
Equation  A-28  into  A-25  giving 


6pK 

°a'  *  °r  2K  (1  +  p?  +  3K 
m  s 


(A-30) 


The  total  axial  stress  is  found  by  substitution  of  A-29  and  A-30 
into  Equation  A-2; 


2Kjj^(l  +  p)  +  6uKg 

2K  (1  +  p)  +  3K^ 
in  s 


(A-31) 


Finally,  the  volume  change  is  obtained  by  plugging  Equation  A-29 
into  A-11; 


AV 


„  2(1  p) 

2K  (1  +  p)  +  3K 
m  s 


(A-32) 


Note  that  accofding  to  Equation  A-23,  the  isotropic  strain, 
e^,  equals  half  the  volume  change,  AV. 

Equations  A-28  through  A-32  are  included  in  the  summary 
of  Table  A-1.  They  are  also  compared  to  the  equations  for  true 
isotropic  loading  derived  from  the  relationships  given  by  Blouin 
and  Kim  (1984). 

PARTIALLY  COUPLED  MODEL 


In  the  partially  coupled  model  an  additional  component  of 
skeleton  volume  strain  due  to  compression  of  the  solid  grains  by 
the  pore  pressure  is  considered.  This  alters  the  last  two 
compatibility  equations  (Equations  A~20  and  A-25)  used  in  the 
derivation  of  the  decoupled  models.  In  addition  to  the  skeleton 
volume  change  due  to  the  mean  effective  stress  given  by  Equation 
A-20,  there  is  a  volume  change  component  due  to  the  pore  pressure 
applied  to  the  solid  grains.  This  component,  AV^,  is  given 
by 


AV 


u 


(A-33) 


for  an  original  skeleton  volume  of  unity.  Thus,  the  skeleton 
component  of  volume  strain  equals  the  volume  strain  in  the 
individual  solid  grains  given  by  Equation  A-7.  In  other  words, 
the  total  skeleton  volume  contracts  in  proportion  to  the 
contraction  in  the  individual  grains.  Assuming  an  isotropic 
skeleton,  the  total  volume  strain  due  to  the  pore  pressure 
loading  will  equal  the  volume  strain  in  the  individual  grains. 
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The  overall  volume  strain  is  obtained  by  combining  Equation  A-33 
with  A-20; 


(A-34) 


The  compatibility  equation  for  zero  axial  strain  is  also 
affected  by  the  skeleton  strain  due  to  the  pore  pressure.  A 
third  term  must  be  added  to  Equation  A-24  to  account  for  the 
additional  axial  strain  in  the  skeleton  from  the  pore  pressure. 
Since  the  volume  strain  is  isotropic,  the  axial  strain  will  equal 
one  third  of  the  volume  strain  from  Equation  A-33,  or 


u 

^a=  3Kg 


Combining  equations  A-35  and  A-24  gives 


(A-35) 


(A-36) 


Setting  the  axial  strain  to  zero  and  putting  Eg  in  terms  of 
K  according  to  Equation  A-17  gives  the  effective  axial 

O 

stress  as 

K^d  -  2’j) 

°a’  "  “  ~~ — K -  (A-37) 

g 
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Thus,  the  five  compatibility  equations  for  the  partially 
coupled  isotropic  loading  in  plane  strain  are 


o  =  o  •  +  u 
r  r 


(A-1) 


(A-2) 


AV= 

m 


(A-11) 


2o  '  +  a  ' 

AV  =  a 

3K^  K 

s  g 


(A-34) 


K^(l  -  2u) 
=  2ua^'  -  u  ^ — - - 


(A-37) 


Solution  of  the  compatibility  conditions  proceeds  as 
follows.  Equating  A-34  and  A-11  gives  the  pore  pressure  as  a 
function  of  the  effective  stresses 


u  (2o^’  +  a^')  3K  (K  K  )  (A-38) 

s  g  m 

Substitution  of  Equation  A-37  into  A-38  gives  the  pore  pressure 
in  terms  of  the  radial  effective  stress  as 
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u  =  2a^' 


(A-39) 


uric 

S 


T^^^nr-rtngr 


Finally,  substitution  of  Equation  A-39  into  A-1  yields  the 
isotropic  effective  stress  in  terms  of  the  applied  total  stress  as 


3K  -  2(1  +  u)K 
_  I  _  _  s  q  ms 

r  “  r  2(1  +  U)  K„(K  -  K  )  +  3K  K 

mg  s  s  9 


(A-40) 


Further  substitution  of  Equation  A-40  into  Equation  A-1  gives  the 
pore  pressure  as  a  function  of  the  applied  total  stress 


2(1  +  y)  K  K 

*  2(1  +  y)K„(K„  -  K^)  +  3K^K„ 

mg  s  s  g 


(A-4I) 


The  effective  axial  stress  is  obtained  by  substituting  Equations 
A-40  and  A-41  into  compatibility  condition  A-37  giving 


6yK  K  -2(l+u)KK 
,  _  ^  s  q  _  ms 

°a  2(1  +  u)K^(K  -  K^)  +  3K 

mg  s  s  g 


(A-42) 


Total  axial  stress  is  obtained  by  substitution  of  A-41  and  A-42 
into  Equation  A-2, 
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(A-43) 


6uK  K  2(1  +  vi)K  (K  -  K  ) 
_  =  _  s  q  m  q  s' 

°a  °r  2(1  +  uiK  (K  -  K  )  +  3K  K 

mg  s  ~ 


Lastly r  the  volume  change  per  unit  volume  is  found  by  substitution 
of  A~41  into  A-11, 


AV 


2(1  +  u)K^ 

2 (1  +  y)K  (K  -  K  )  +  3K  K 
mg  s  s  9 


(A-44) 


Recall  that  from  Equation  A-23,  the  isotropic  strain,  e^, 
simply  equals  half  the  volume  strain  given  by  Equation  A-44. 

The  stress,  pore  pressure  and  volume  strain  relationships 
from  Equations  A-40  through  A-44  for  the  partially  coupled  model 
are  summarized  in  Table  A-1,  They  are  compared  to  the  Equations 
for  a  true  isotropic  loading  derived  using  the  partially  coupled 
model  by  Blouin  and  Kim  (1984). 


FOLLY  COUPLED  MODEL 


The  fully  coupled  model  utilizes  all  the  assumptions  of 
the  partially  coupled  model  and  in  addition  takes  into  account 
the  volume  change  in  the  solid  grains  due  to  application  of  the 
effective  stresses  in  the  skeleton.  Of  the  five  compatibility 
equations  listed  in  the  previous  subsection  (Equations  A-1,  A-2, 
A-11,  A-34  and  A-37)  only  Equation  A-11  is  modified  by  this  added 
compatibility  condition.  The  other  volume  change  equation,  A-34, 
is  not  modified  because  the  definition  of  the  skeleton  bulk 
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modulus*  Kg,  automatically  takes  account  of  the  solid  grain 
strains  resulting  from  the  effective  stresses. 

The  effective  stresses  are  defined  in  terms  of  the  total 
cross  sectional  area  of  a  given  material  element.  Because  only  a 
portion  of  the  cross  sectional  area  is  made  up  of  solid  grains* 
the  actual  stress  in  the  solid  grains*  termed  the  intragranular 
stress*  is  higher  than  the  effective  stress.  The  intragranular 
stresses  are  obtained  by  dividing  the  effective  stresses  by  (1  -  n)* 
the  fraction  of  the  total  volume  occupied  by  the  solid  grains. 

Thus*  the  intragranular  isotropic  stresses*  Oj-j^'*  and  the 
intragranular  axial  stress*  a^^**  are  given  by* 


and 


(A-45) 


The  mean  intragranular  stress*  o  ,  is  the  average  of  the 
sum  of  the  two  components  of  isotropic  stress  and  the  axial 
stress*  given  by 


2o  '  +  a  ’ 

=  ■  3(--  n? 

The  volume  strain  in  a  unit  volume  of  saturated  porous  material 
due  to  the  intragranular  stress  is  obtained  by  dividing  the  mean 
stress  by  the  bulk  modulus  of  the  solid  grains  and  multiplying  by 
the  volume  of  the  solid  grains  in  the  unit  element*  (1  -  n). 
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ThuSf  the  volume  strain  due  to  the  intragranular  stress,  AV^, 
is  given  by 


AV, 


3Ti“  n)Kg 


(A-47) 


or  simply 


AV, 


+  o. 


(A-48) 


The  total  volume  change  in  the  saturated  unit  element  is  obtained 
by  combining  the  volume  change  due  to  the  intragranular  stress 
given  by  Equations  A-48  with  the  volume  change  due  to  the  pore 
pressure  in  the  pore  fluid  and  on  the  solid  grains  given  by 
Equation  A-11.  The  total  volume  change  is  thus  given  by 


„  2o^'  +  o  • 

“  ir  *  ■■ 

m  g 


(A-49) 


The  five  compatibility  equations  for  the  fully  coupled 
plane  strain  loading  are  given  by 


a  =  a  '  +  u 
r  r 


(A-1) 


a,  =  a  '  +  u 
a  a 


(A-2) 
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LV  = 


(A-49) 


2o  '  +  o  ' 
u  .  r  a 

K  — Sir - 

m  g 


--jK  —  ^  ir 

s  g 


K^(l  -  2y) 


“a'  =  •  “  K 


Equating  Equations  A-49  and  A-34  gives  the 
in  terms  of  the  effective  stresses  as 


2o^'  +  o,'  K„(K„  -  K^) 

r  a  m  g  s 

3  K^(K„  -  K  ) 

s  g  m 


For  simplicity  Equation  A-50  is  reexpressed  as 


u  = 


^'’r'  *  “a' 


Where  R|^  is  the  modulus  ratio  given  by 


p  —  ^  9  ^ 

^  K.(K  -  K  ) 


m 


Substitution  of  Equation  A-37  into  A-51 
effective  isotropic  stress,  o  j- ' »  as  a  function  of 
pressure , 


(A-34) 

{A-37) 

pore  pressure 


(A-50) 


(A-51) 


(A-52) 

the 

pore 
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3K^  +  (1  -  2  u )K^R^ 
2(1  +  m)K^R^ 


(A-53) 


Substitution  of  A-53  into  Equation  A-1  gives  the  pore 
pressure  in  terms  of  the  applied  total  isotropic  stress  as 


(A-54) 


where  the  denominator  is  given  by 

Df  =  \  (2(1  +  u)Kg  +  (1  -  2u)Kg 

Substitution  of  Equation  A-54  into  a-53  gives  the  isotropic 
effective  stress  as 


) 


+  3K 


(A-55) 


Cr’  = 

3K^  +  (1  -  2  U  )KgRj^ 

ective  stress  is  obtained  by  substi 

into  A-39 

to  give 

6uK^  -  2(1  -  2u)KgI^ 

®a 

(A-56) 


(A-57) 


The  total  axial  stress  is  obtained  by  substituting  Equations  A-54 
and  A-57  into  A-2  to  give 
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2Rj^  (d  +  u)K^  -  (1  -  2u)Kg 

+  6uK^ 

“a  "r  Dj 

(A-58) 


Finally,  the  volume  strain  in  the  unit  element  is  found  by 
substitution  of  Equations  A-54,  A-56  and  A-57  into  Equation  A-34 
giving 


(A-59) 


The  fully  coupled  plane  strain  relationships  between 
effective  stresses,  pore  pressure,  volume  strain  and  the  applied 
total  isotropic  stress  are  summarized  in  Table  A-1.  They  are 
compared  to  corresponding  relationships  for  a  true  isotropic 
loading  derived  from  the  report  by  Blouin  and  Kim  (1984). 


NUMERICAL  EXAMPLE 

A  numerical  example  illustrating  the  differences  between 
the  various  models  of  Table  A-1  is  listed  in  Table  A-2.  Material 
properties  of  a  typical  saturated  sandstone  having  a  porosity  of 
20%  were  assumed;  these  properties  are  listed  in  Table  A-2.  The 
axial  symmetric  loading  in  plane  strain  is  compared  to  the 
isotropic  loading  for  each  of  the  three  models.  There  are 
significant  differences  between  the  plane  strain  and  isotropic 
loadings  as  well  as  significant  differences  between  the  various 
models . 
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For  each  model,  the  volume  change  computed  for  the 
isotropic  loading  is  significantly  (14.0%  to  17.7%)  higher  than 
the  volume  change  under  the  plane  strain  loading.  The  pore 
pressures  developed  under  the  isotropic  loading  are  also 
correspondingly  higher  than  those  from  the  plane  strain  loading. 
The  effective  isotropic  stresses,  however,  are  about  10%  to  14% 
less  than  the  corresponding  effective  stresses  in  the  plane 
strain  loading. 

There  are  also  significant  differences  between  the 
various  models.  For  the  plane  strain  loading,  the  volume  change 
for  the  fully  coupled  model  is  22.3%  larger  than  that  for  the 
decoupled  model  and  8.2%  larger  than  that  for  the  partially 
coupled  model.  A  similar  trend,  with  even  larger  differences, 
holds  for  the  isotropic  loading.  In  general,  the  pore  pressures 
computed  using  the  partially  coupled  model  are  higher  than  those 
from  the  fully  coupled  model,  while  the  pressures  from  the 
decoupled  model  are  lower  than  those  from  the  fully  coupled 
model.  The  effective  stresses  from  both  the  decoupled  and 
partially  coupled  models  are  less  than  those  from  the  fully 
coupled  model,  though  the  decoupled  effective  stresses  are  in 
closer  agreement  than  are  the  partially  coupled  effective 
stresses. 
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Table  A.l.  Comparison  of  axisymmetric  and  isotropic  loadings 


TeUsle  A. 2.  Numerical  comparison  between  undrained  axial 

symmetric  ^Tjlame  strain)  and  isotropic  loadings. 


Material:  Saturated  Sandstone 

Properties:  Bulk  Modulus  of  Water 

Bulk  Modulus  of  Solid  Grains 

Bulk  Modulus  of  Skeleton 

Porosity 

Poisson's  Ratio 

Mixture  Modulus  (Eq.  A-13) 


-  0.29  E6  psi 
Kg  *  5  E6  psi 
Kg  ■  1.5  E6  psi 
n  -  0.2 

li  -  0.2 

K  »  1.177  E6  psi 

IB 


MODEL 

AXIAL  SYMMETRIC  LOADING 

IN  PLANE  STRAIN 

ISOTROPIC  LOADING 

Decoupled  Model 

Pore  Pressure 

u  •  0.386  0^ 

u  *  0.440 

Effective  Stress 

Oj.’  =  0.614  0 

0  '  *  0.560  a 

Total  Axial  Stress 

Effective  Axial  Stress 

Volume  Change 

0  »  0.631  0 

0.'  -  0.246  0^ 

av  X.  3.277  E-7  0^ 

. 

av  .  3.736  E-7  Oy 

Partially  Coupled  Model 

Pore  Pressure 

u  «  0.436  0 

u  .  0.506  a 

Effective  Stress 

r 

a*  »  0.564  0^ 

r 

»  0.494  0 

Total  Axial  Stress 

Effective  Axial  Stress 

Volume  Change 

r  r 

o_  .  0.583  a_ 

A  ^ 

0^'  «  0.147  0^ 

A  £ 

av  =  3,705  E-7  0^ 

r  r 

av  »  4.303  E-7  0 

r 

Fully  Coupled  Model 

Pore  Pressure 

u  »  0.355 

u  =  0.418  0 

Effective  Stress 

r 

0  '  «  0.645 

r 

0  '  *  0.582  0 

Total  Axial  Stress 

Effective  Axial  Stress 

Volume  Change 

r  r 

0,  =  0.549  o_ 

3  t 

a'  =  0.194  o_ 

A  Jm 

av  »  4.007  E-7  0^ 

r  r 

av  «  4.716  E-7  0 

r 

A24 


Total  Si 
Effectl-' 


APPENDIX  B 


FLOW  DUE  TO  VOLUME  CHANGE  IN  SATURATED  POROUS  MATERIALS 


A  fundamental  aspect  of  flow  in  saturated  porous  media  is 
the  relationship  between  pore  water  flow  and  volume  change.  This 
relationship  is  the  basis  for  accurate  flow  calculations  in 
multi-phase  computational  models  and  is  used  in  CONSL  as  one 
method  of  calculating  flow  (see  Section  5).  For  porous  fully 
saturated  materials  which  obey  the  effective  stress  law,  i.e.  the 
total  stress  on  an  element  of  material  equals  the  pressure  in  the 
pore  water  throughout  the  material  plus  the  intergranular  or 
effective  stress  in  the  material  skeleton,  the  flow  of  pore  water 
into  or  out  of  any  element  of  material  due  to  an  increase  or 
decrease  of  stress  and/or  pore  water  pressure  in  the  element 
equals 


—  the  change  in  volume  of  the  element,  as  measured 
by  change  in  skeleton  volume, 

—  plus  pore  water  flow  into  or  out  of  the  element 

due  to  compression  or  expansion  of  the  pore  water  and 
solid  grains  within  the  element. 

The  total  flow  of  pore  water,  AV^,  is  given  by 

AV^  =  AVg  -  AVg  -  AV^  (B-1) 


where  AVg  is  the  change  in  skeleton  volume  of  the  element. 


B1 


AVg  is  the  change  in  volume  of  the  solid  grains  in  the  element 
due  to  an  increase  or  decrease  in  pore  pressure,  and  is  the 

change  in  volume  of  the  pore  water  in  the  element  due  to  an 
increase  or  decrease  in  pore  pressure.  The  signs  indicate  that 
flow  due  to  change  in  skeleton  volume  is  opposite  to  flow  due  to 
volume  changes  within  the  skeleton.  For  example,  an  increase  in 
total  stress  results  in  compression  of  the  skeleton  with  flow  of 
pore  water  out  of  the  matrix,  while  the  corresponding  increase  in 
pore  pressure  within  the  element  compresses  the  pore  water  and 
solid  grains  requiring  a  component  of  flow  into  the  element. 

Pore  water  flow,  as  computed  from  Equation  B-1,  will 
depend  on  the  coupling  relationships  assumed  in  the  material 
models.  Blouin  and  Kim  (1984)  describe  three  coupling 
relationships  between  the  material  skeleton  and  pore  water  having 
varying  degrees  of  accuracy  depending  on  the  level  of 
sophistication  desired.  These  are  also  described  in  Appendix  A. 
Flow  equations  for  each  of  the  three  models  are  developed  in  the 
following  sections. 


DECOUPLED  MODEL 


The  decoupled  model  is  the  simplest,  though  least 
accurate,  of  the  three  models. 

The  decoupled  model  includes  compressible  pore  water  and 
solid  grains,  but  the  skeleton  and  pore  water  are  assumed  to  act 
in  parallel,  with  no  coupling  between  the  two.  The  skeleton 
volume  strain  is  given  by 


AV 


sd 


Ao,^' 

m 

K 


82 


(B-2) 


where  AQjjj'  is  the  change  in  mean  effective  stress  and  Kg 
the  bulk  modulus  of  the  material  skeleton.  The  change  in  volume 
of  the  solid  grains  within  the  soil  element  due  to  compression  by 
the  pore  water  is  given  by 


(1  -  n) 


(B-3) 


where  n  is  the  porosity,  Au  the  change  in  pore  pressure,  and 
Kg  the  bulk  modulus  of  the  solid  grains.  The  term  (1  -  n) 
represents  the  volume  fraction  of  solid  grains.  The  change  in 
volume  of  the  pore  water  in  the  element  due  to  the  change  in  pore 
water  pressure  is  given  by 


w 


(B-4) 


where  K^  is  the  bulk  modulus  of  the  pore  water.  Substitution 
of  Equations  B-2,  B-3  and  B-4  into  Equation  B-1  gives  the  flow 
volume  for  the  decoupled  model  as 


Ao_’ 

AVf^  =  -  AU 


(B-5) 


This  can  be  rewritten  as 
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(B-e; 


Ao 


AV 
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fd 


-  iiu 


(K  +  nK  -  nK  \ 
w _ w _ a  ) 

K  K  / 

q  w  • 


The  term  in  the  parentheses  is  the  compressibility  of  the  solid 
grain/pore  water  mixture  derived  in  Appendix  A  and  given  in 
Equation  A-13.  Equation  B-6  can  be  reexpressed  as 


AV 


fd 


Aa_' 

m 

K_ 


(B-7) 


where  is  the  bulk  modulus  of  the  solid  grain/pore  water 
mixture  with 


K  K 

K  -  _ g__w _ 

m  K  +  n(K  -  K  ) 
w  g  w 


(B-8) 


Thus,  the  flow  in  the  decoupled  model  simply  equals  the  change  in 
volume  of  the  material  skeleton  due  to  a  change  in  effective 
stress,  less  the  change  in  volume  of  the  solid  grain/pore  water 
mixture  due  to  a  change  in  pore  pressure. 


PARTIALLY  COUPLED  MODEL 

In  the  partially  coupled  model  the  change  in  skeleton 
volume  due  to  compression  of  the  grains  by  the  pore  water 
pressure  is  included  in  addition  to  the  relationships  in  the 
decoupled  formulation.  The  only  change  in  the  volume  change 
compatibility  equations  in  the  partially  coupled  model  is  in  the 
skeleton  volume  change  equation.  In  the  partially  coupled  model 
the  skeleton  volume  is  a  function  of  both  the  applied  effective 
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stress  and  the  compression  of  the  grains  making  up  the  skeleton 
by  the  applied  pore  pressure.  The  partially  coupled  skeleton 
volume  strain  is  given  by 


(B-9) 


The  last  term  in  Equation  B-9  is  the  skeleton  volume  strain  due 
to  the  applied  pore  pressure  increment  Au.  As  explained  in 
Appendix  A,  the  total  skeleton  volume  strain  equals  the  volume 
strain  in  the  individual  particles  (assuming  isotropic  material). 
The  remaining  volume  change  compatibility  equations  are  identical 
to  those  of  the  decoupled  model  given  in  Equations  B-3  and  B-4. 
Note  that  the  volume  change  of  the  solid  grains  given  by  B-3  is 
included  in  the  decoupled  model,  even  though  this  influence  is 
not  included  in  the  volume  change  of  the  skeleton. 


Substitution  of  Equations  B-3,  B-4  and  B-9  into  Equation 
B-1  gives  the  partially  coupled  flow  volume  as 


(B-10) 


In  terms  of  the  bulk  modulus  of  the  solid  grain/pore  water 
mixture.  Equation  B-10  is  expressed  as 


Au 

K 


m 


(B-11) 
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where  the  change  in  skeleton  volume  AV__  is  given  by  Equation 
B-9. 


FULLY  COUPLED  MODEL 


The  fully  coupled  model  described  by  Blouin  and  Kim 
(1984)  and  in  Appendix  A  is  the  most  sophisticated  of  the  three 
multi-phase  models.  The  fully  coupled  model  incorporates  all  the 
features  of  the  partially  coupled  model  plus  it  accounts  for 
volume  change  in  the  solid  grains  resulting  from  the  effective 
stresses  in  the  grains.  This  additional  compatibility  condition 
alters  only  Equation  B-3  of  the  partially  coupled  compatibility 
equations.  In  the  fully  coupled  model  the  volume  of  water 
flowing  in  or  out  of  the  porous  element  due  to  compression  of  the 
solid  grains  is  made  up  of  a  component  due  to  compression  by  the 
pore  pressure  and  a  component  due  to  compression  of  the  solid 
grains  by  the  effective  stresses  in  the  skeleton.  The  overall 
volume  change  in  the  grains  is  given  by 


(1  -  n) 


(B-12) 


where  the  first  term  is  the  volume  compression  by  the  pore  water 
as  given  in  Equation  B-3.  The  last  term  is  the  compressive 
volume  strain  in  the  solid  grains  due  to  the  mean  effective 
stress  increment  The  actual  stress  in  the  grains  is 

termed  the  intragranular  stress.  This  is  higher  than  the 
intergranular  or  effective  stress  because  the  effective  stress  is 
an  average  value.  In  actuality  the  total  effective  loading  is 
carried  only  in  the  grain  structure,  not  in  the  intervening  pore 
spaces.  Thus,  the  mean  intragranular  stress  is  given  by 
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The  volume  change  due  to  the  intragranular  stress,  AVg^  is 


Ao 


m  1 


(1  -  n) 


Substitution  of  Equation  B-13  into  B-14  gives 


{B-14) 


m 


(B-15) 


which  is  the  last  term  of  Equation  B-12.  Note  that  the 
intragranular  stress  does  not  enter  into  the  volume  change 
relationships  for  the  soil  skeleton  given  by  Equation  B-9.  This 
is  because  the  intragranular  stresses  are  automatically  accounted 
for  in  the  skeleton  bulk  modulus.  Kg. 


Combining  the  fully  coupled  volume  change  components  from 
Equations  B-4,  B-9  and  B-12  into  Equation  B-1  gives  the  volume 
of  water  flowing  into  or  out  of  a  fully  coupled  element  as 


(B-16) 


In  terms  of  the  mixture  bulk  modulus,  flow  is  given  as 
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where  AV^^  is  the  volume  change  of  the  skeleton  given  by 
Equation  B-9. 


In  numerical  two  phase  calculations  output  is  often  in 
terms  of  volume  change  and  pore  pressure.  Flow  volume  can  be 
calculated  directly  from  these  parameters  by  a  slight 
modification  of  Equation  B-17.  Manipulation  of  Equation  B-9 
gives  the  change  in  mean  effective  stress  as 


La'  =  -  Au  ^  (B-IS) 

m  s  s  Kg 

Substitution  of  Equation  B-18  into  B-17  yields 


CHECK  CASE  -  UNDRAINED  LOADING 


As  a  check  on  the  fully  coupled  flow  Equation  B-17,  the 
undrained  effective  stress  and  pore  pressure  relationships  derived 
for  the  fully  coupled  model  in  Appendix  A  can  be  used  to  certify 
that  there  is  no  flow  into  or  out  of  a  saturated  element  during 
an  undrained  loading.  Substitution  of  Equation  B-9  into  Equation 
B-17  gives 
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For  the  plane  strain  hydrostatic  loading  condition,  the  mean 
effective  stress  is  given  by 


(B-21) 


where  Oj- '  is  the  radial  hydrostatic  effective  stress  and 
Og '  is  the  axial  effective  stress  in  which  direction  no  strain 
is  allowed.  Substitution  of  the  effective  stresses  from  Equations 
A-56  and  A-57  into  B-21  gives  the  mean  effective  stress 


m 


•  =  o 


2Kg(l  +  w) 

r  d" 


(B-22) 


in  terms  of  the  total  stress,  Poisson's  ratio  for  the 

soil  skeleton  is  y  and  the  denominator  is  given  by 
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Substitution  of  the  mean  effective  stress  from  Equation  B-22  and 
the  pore  pressure  from  Equation  A-54  into  B-20  gives 
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Substitution  of  Equation  B-25  into  B-24  and  further  manipulation 
gives 


”  K  K  (K  -  K  i 
m  s  g  m 


(B-26) 


Thus,  the  computed  flow  is  zero,  which  satisfies  undrained 
loading  compatibility  conditions. 


NOTE  OF  CAUTION 


In  applying  the  flow  equations  based  on  the  skeleton 
volume  change,  care  must  be  taken  not  to  violate  the 
compatibility  conditions  implicit  in  the  flow  equations.  In 
making  finite  element  calculations,  it  is  possible  to  specify 
boundary  and  pore  pressure  conditions  which  violate  the 
compatibility  conditions  used  to  derive  the  equations.  This 
occurs  when  boundary  conditions  are  set  so  that  the  deformations 
of  the  individual  elements  of  the  calculational  grid  do  not 
correspond  to  the  deformation  of  the  material  skeleton.  For 
instance,  if  grid  boundaries  are  fixed  and  pore  pressure  is 
applied  to  the  saturated  porous  material,  the  material  skeleton 
will  contract  due  to  grain  compression  by  the  pore  water.  Since 
the  grid  boundaries  are  fixed,  however,  the  element  boundaries 
will  not  contract  with  the  material  skeleton.  In  this  case  large 
errors  in  any  flow  calculations  across  the  grid  boundaries  are 
introduced. 
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APPENDIX  C 


FINITE  ELEMENT  FORMULATION  OF  CONSOLIDATION 


It  is  a  common  assumption  in  soil  mechanics  that  the 
solid  grains  and  pore  water  making  up  a  saturated  soil  are 
incompressible.  This  assumption  is  a  reasonably  good 
approximation  for  soils  whose  skeleton  moduli  are  much  lower  than 
those  of  the  solid  grains  and  pore  water.  For  most  rocks, 
however,  the  skeleton  modulus  is  on  the  same  order,  or 
significantly  greater  than  that  of  the  pore  water.  Thus,  the 
compressibility  of  the  solid  grains  and  pore  water  plays  an 
important  role  in  the  consolidation  process  of  saturated  porous 
rocks. 


Presented  in  this  appendix  is  the  finite  element 
formulation  of  consolidation  including  the  compressibility  of 
solid  grains  and  pore  water. 


NOTATION 


Note  that  positive  signs  have  been  used  for  elongation 
and  tension.  A  comma  denotes  differentiation  with  respect  to  the 
subsequent  indices  and  the  superposed  dot  denotes  time  rate. 

{u}  :  solid  phase  displacement 

(U)  :  fluid  phase  displacement 

{a}  :  total  stress 

{a'}  :  effective  stress 

ir  :  fluid  pressure 


Cl 


{ e  } 

(u}e 

(и) 
{n} 
{T} 

yv 

Q 

{b} 

[к] 

[oeP] 

(1} 


C 

s 

a 

P 

Pf 

^ij 
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[C] 

[E] 

[H] 

(f) 

(r) 

(Q) 


solid  phase  strain 

solid  phase  volumetric  strain 

fluid  phase  volumetric  strain 

nodal  solid  phase  displacement  vector  at  the 

element  degrees  of  freedom 

nodal  fluid  pressure  at  the  element  degrees  of 
freedom 

nodal  solid  phase  displacement  vector  at  the 

structural  degrees  of  freedom 

nodal  fluid  pressure  at  structural  degrees  of 

freedom 

applied  boundary  traction 

specified  boundary  flux 

body  force  vector 

permeability  matrix 

elasto  plastic  stress-strain  matrix 

unit  vector  (1)^=  <11I  00  0> 

porosity 

fluid  compressibility 
compressibility  of  solid  grains 
compressibility  of  soil-water  mixture 
bulk  mass  density  of  mixture 
fluid  mass  density 
Kronecker's  delta 
tangent  stiffness  matrix 

coupling  matrix  between  solid  skeleton  and  pore 
f  luid 

matrix  of  compressibility  of  pore  fluid 
dissipation  resistance  matrix 
vector  of  nodal  forces 
internal  resisting  force  vector 
equivalent  flow  vector 
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FIELD  EQUATIONS 


Principle  of  effective  stress 


°ij  “  ‘^ij'  °ij  ^ 


(C-1) 


Constitutive  law  for  solid  skeleton 


{do'}  =  [D®P]  (del  - 


-f  {1}  dTT 


(C-2) 


Continuity  Equations 

The  coupled  continuity  equation  of  flow  has  been  derived 
by  Kim,  1982a  and  is  given  by 


(1  -  n)de^  +  ndEp  -  adir  -  (1  -•  n)Cg'  dp  =0  (C-3) 


where 


a=nC  +(l-n)C 
w  s 
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C  '  =  - = - 

s  (1  -  n) 
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Substitution  of  Equation  C-5  into  Equation  C-3  gives 


n(de  -  de  )  =  adTT  +  C  dp'  -  de  (C-6) 

r  V  s  V 


Generalized  Darcy's  Flow  Law 

Neglecting  inertia  term. 


w^  =  n(U^  >  u^)  =  (Tr,j  +  p^bj)  {C-7) 

FINITE  ELEMENT  DISCRETIZATION 

Two  global  equilibrium  equations  are  derived,  first  in 
terms  of  field  variables  and  then  discretized  using  the  nodal 
variables. 

The  total  stresses  are  in  equilibrium  with  the  applied 
boundary  tractions.  Taking  the  solid  skeleton  movement  as  the 
virtual  displacement,  the  internal  and  external  virtual  work  must 
be  equal.  At  certain  time  t. 
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dv 
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and 


5W^  =  5W_  (C-10) 

The  internal  fluid  movements  relative  to  the  solid 
skeleton  are  compatible  with  the  specified  boundary  flux.  Taking 
the  fluid  pressure  field  at  a  certain  time  t  as  the 
complementary  virtual  stresses,  the  internal  and  external 
complementary  virtual  work  must  be  equal. 


and 


5Wj*  *  5Wg* 


(C-14) 


Discretization  of  field  variables  into  nodal  variables 


{u}  =  [n]  {u}g 
ic)  =  [B]  in)^ 
IT  =  <G> 

{tt,  J  =  [A]  {i}g 


FINITE  INCREMENTAL  RELATIONSHIPS 


Stress  vector  at  time  step  n  can  be  expressed  as 


(a  }  =  (a  '  ,  }  +  {Ao„'}  +  {1}  IT  (C-16) 

n  n-1  n  n 


Substitution  of  Equation  C-2  into  Equation  C-16  yields 
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(0„)  -  [D®P]  tic„)  +  ^(1)  -  -j 


:„>  +  (fl)  -  T  'l’)'n 


Substituting  Equation  C-17  into  Equation  C— 8  and  replacing  the 
field  variables  in  Equations  C-8  and  C-9  by  the  nodal  values 
using  Equation  C-15  gives. 


=  {6u} 
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[N]^p  (b)  dv  +  zl  [N]^  {T}  ds 
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From  Equation  C-10 


[K^]  taO^)  +  Cc]  {i„)  =  (F^)  - 
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where 
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The  right  hand  side  of  Equation  C-6  can  be  expressed  as 


ad-rr  +  dp'  -  de^  =  -  -|—  {l}*^  [D®P 

-  i  {l}"^  [D®P]^{de} 


Substituting  Equation  C-22  into  Equation  C-12  and 
replacing  the  field  variables  in  Equations  C-12  and  C-13  by  the 
nodal  values  using  the  Equation  C-15, 
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From  Equation  C* 
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In  deriving  Equation  C~25,  the  symmetry  of  [D®P]  has  been 
assumed.  Assuming  the  linear  variation  of  skeleton  displacements 
and  fluid  pressures  between  t  and  t  At  and  setting  the 
equilibrium  equations  at  time  t  +  l/2Atr 


[C]*^ 


{W}  =  {Q} 
n  n 


(C-27) 


where 


f5„)  =  -  ^  ((Q„)  *  ^  [H]  -  [E])  (C-28) 

Equations  C-20  and  C-27  can  be  solved  simultaneously  by  the 
incremental  step  by  step  procedure.  These  incremental  equations 
have  been  implemented  in  the  one-dimensional  axisymmetric  finite 
element  code  "CONSL. " 


CIO 


onnonnono 


APPENDIX  D 


CONSL  CODE  LISTING 


FTN7X,E 

$FILES(0,5) 

$EMA  /WARE A/ 

PROGRAM  CONSL 


REVISION: 

REVISION: 

REVISION: 

REVISION: 

REVISION: 


August  30,  1984 
September  1,  1984 
September  2,  1984 
September  5,  1984 
September  6,  1984 


INTEGER  FILE1(6) ,FILE2(6) 

COMMON/WAREA/A( 32767 ) ,M ( 14000 ) 

COMMON/GINPT/TITLE( 20) , RI , RO,NUMEL, NKF , NEF , NTF , IPLANE 
COMMON/GPARM/NLNR, NPK , NPE , NPH ,NUMNP , NEO , NEQH , AP , NVF , NDF 
COMMON/BPRES/NMOD,NOD( 2) ,SPR(2) 

COMMON/BFORC/NMOS , NOS ( 2 ) , TNP ( 2 ) , NMOU , NOU ( 2 ) , SPU ( 2 ) 
COMMON/PRINT/TPRNT ( 20 ) , NPRNT 
COMMON/INIST/SRI,STI ,PI 


2021 


I 

I 


C 


c 

c 


c 

c 


WRITE(l)  '  TYPE  NAME  OF  INPUT  FILE 
READ( 1,2021)  FILEl 

WRITE(l)  '  TYPE  NAME  OF  OUTPUT  FILE 
READ(1,2021)  FrLE2 
FORMAT (6A2) 

OPEN(6,FILE=FILEl,STATUS='OLD' ) 
OPEN(7,FILE=FILE2,STATUS='NEW' ) 
OPEN(8,FILE='CONOUT: :-19' ) 
OPEN(9,FILE='FLOW: :-19'  ) 
OPEN(5,FILE='XFLOW: :-19' ) 

NTOTL  =  32767 
ITOTL  =  14000 

CALL  DREAD 

11  =  1 

12  =  I1+2*NUMNP 

13  =  I2+NUMNP 

14  =  I3+NUMNP 

15  =  I4+2*NUMEL 

16  =  I5+2*NUMEL 

17  =  I6+2*NUMEL 
N1  =  1 

N2  =  Nl+NUMNP 
CALL  NODEG(A(Nl ) ) 


I 


I 
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CALL  ELEMT (M(I1),M(I2),M(I3),M{I4),M(I5),M(I6), NUMEL , 

*  NUMNP,NE0,NEQH,A(N1 ) ) 

C 

18  =  I7+NE0 

19  =  I8+2*NE0H 

IF( I9.GT.IT0TL)  GO  TO  200 
N3  *  N2+4*NEO 
N4  -  N3 

IF(NLNR.NE.O)  N4  =  N3+4*NEO 
N5  =  N4+2*NEOH 
N6  =  N5 

IF(NLNR.GT. 1. AND.NLNR.lt. 6)  N6  =  N5+2*NEOH 

N7  =  N6+NEOH 

N7A=N7 

IF(NVF.EO. I.AND.PI.NE.0.0)  N7A=N7+NEQH 

N8  =  N7A+NEOH 

N9  =  N8+NEQ 

NIO  =  N9+NEO 

Nil  =  NIO+NEQ 

N12  =  Nll+NEO 

N12A  =  N12 

IF(NMOD.EO.l.OR.NMOU.EO.l)  N12A  =  N12+NEO 
C 

N13  =  N12A+2*2*NUMEL 
N14  =  N13 

IF(NKF.EO.l)  N14  =  N13+2*NUMEL 
IF(N14.GT.NTOTL)  GO  TO  200 
C 

NEOH2  =  2*NE0H 
C 

CALL  KSTAR(A(N1) ,A(N2) , A ( N3 ) ,M ( 17 ) , A( N4 ) , A ( N5 ) , M ( I 8 ) , 

*  A(N6),A(N11) ,M(I4),M(I5),M(I6) ,A(N12A), 

*  A(N13) ,NUMEL,NUMNP,NEO,NEOH,NEOH2,A(N12) ,A(N7) ) 
C 

CALL  SOLVE(M(I3) ,A(N1) ,A(N2) ,A{N3) ,M(I7) ,A(N4) ,A(N5) , 

*  M(I8) ,A(N6) ,A(N7A) , A( N8 ) , A ( N9 ) , A ( NlO ) , 

*  A(Nll) ,M(I4) ,M(I5) ,M(I6) ,A(N12A) ,A(N13) , 

*  NUMEL, NUMNP,NEO,NEOH, NEOH2, A( N12 ) ,A(N7) ) 

C 

200  CONTINUE 

WRITE( 7,2000)  ITOTL,I9 

2000  FORMAT ( 'MAX.  INTEGER  DIMENSION  SET  IN  PROGRAM  =  ',15/ 

.  'REQUIRED  INTEGER  DIMENSION  -  =  ',15//) 

C 

WRITE(7,2001)  NT0TL,N14 

2001  FORMAT('MAX.  REAL  DIMENSION  SET  IN  PROGRAM  =  ',15/ 

.  'REQUIRED  REAL  DIMENSION  -  =  ',15//) 

C 

CLOSE(5) 

CLOSE{6) 

CLOSE(7) 

CLOSE(8) 
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CLOSE ( 9 ) 

STOP 

END 


C 

$EMA  /WAREA/ 

BLOCK  DATA 

COMMON /WAREA/A(32767) ,M(14000) 

C0MM0N/GINPTAITLE( 20) ,RI,RO,NUMEL,NKF,NEF,NTF,IPLANE 
COMMON /GPARM/NLNR , NPK , NPE , NPH , NUMNP , NEO , NECH , AP , NVF , NDF 
COMMON/MTPRT/CM , V , POR , PK , ALPHA , SO , STAW , RSG 
COMMON/TIMEV/DT,TEND,NDT,DTT( 10) ,NCL( 10) 
COMMON/INIST/SRI,STI ,PI 

COMMON/BFORC/NMOS , NOS ( 2 ) , TNF (2) , NMOU , NOU ( 2 ) , SPU ( 2 ) 
COMMON/BPRES/NMOD , NOD ( 2 ) , SPR ( 2 ) 

COMMON/BCONC/ JSDL , JDCL , TNFL , BPPL , JSDR , JDCR , TNFR , BPPR , 
#SSDL,SSDR 

COMMON/PRINT/TPRNT(20) ,NPRNT 
END 
C 

SUBROUTINE  DREAD 

COMMON /GINPT/TITLE ( 20 ) , RI , RO,NUMEL, NKF , NEF , NTF , I PLANE 
COMMON/GPARM/NLNR , NPK , NPE , NPH , NUMNP , NEQ , NEQH , AP , NVF , NDF 
COMMON/MTPRT/CM , V , POR , PK , ALPHA , SO , STAW , RSG 
COMMON AIMEV/DT, TEND, NDT,DTT( 10) ,NCL( 10) 
COMMON/INIST/SRI ,STI , PI 

COMMON/BFORC/NMOS , NOS ( 2 ) ,TNF { 2 ) , NMOU , NOU ( 2 ) , SPU ( 2 ) 
COMMON/BPRES/NMOD , NOD ( 2 ) , SPR { 2 ) 

COMMON/BCONC/JSDL, JDCL,TNFL, BPPL, JSDR, JDCR, TNFR, BPPR 
#,SSDL,SSDR 

COMMON/PRINT/TPRNT(20) ,NPRNT 
PROBLEM  SPECIFICATION 

READ(6,1001)  TITLE 

1001  FORMAT (20A4) 

WRITE(7,2001)  TITLE 

2001  FORMAT (/20A4///) 

GLOBAL  INPUT 

READ (6,1002)  RI , RO , AP , NUMEL , NKF , NEF , NTF , I  PLANE , NVF , N DF 

1002  FORMAT (3E10. 0,715) 

WRITE (7,2002)  RI , RO , AP , NUMEL, NKF , NEF , NTF , I PLANE , NVF , NDF 


2002  FORMAT (  'RADIUS  OF  TUNNEL - : - -  ',E10.3/ 

.  'RADIUS  OF  REMOTE  BOUNDARY  -  =  ',E10.3/ 

.  'MESH  GROWTH  FACTOR  -  =  ’,E10.3/ 

.  'NUMBER  OF  ELEMENTS - =  ',15  // 

.  '.EO.O  ELASTIC  SKELETON  '  / 

.  '.EO.l  NONLINEAR  SKELETON  -  =  ',15  // 

.  '.EO.O  SATURATED  FLUID  •  / 

.  '.EO.l  PARTIALLY  SATURATED  FLUID  -  =  ',15  // 

'.EO.O  CONSTANT  TIME  STEPS  '  / 
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.  '.EQ.l  VARIABLE  TIME  STEPS - -  ',15  // 

’.EQ.O  1-D  CYLINDRICAL  SYMMETRY  '  / 

'.EQ.l  1-D  PLANE  STRAIN - =  ',15  // 

.  '.EO.O  FLOW  VOL.  BASED  ON  FIRST  ELEM  '  / 

.  '.EQ.l  FLOW  VOL.  BASED  ON  ALL  ELEM.-  =  ',15  // 

.  '.EO.O  DECOUPLED  STORAGE  EQUATION  '  / 

.  '.EQ.l  FULLY  COUPLED  STORAGE  EQUATION=  ',15  ///) 

C 

NUMNP  =  NUMEL+1 


INPUT  FOR  MATERIAL  PROPERTIES 


READ(6,1003)  CM,V,POR,PK 
1003  FORMAT(E10.3,2F10.0,L10.3) 

WRITE(7,2003)  CM,V,POR,PK 

2003  FORMAT ( 'SKELETON  CONSTRAINED  MODULUS  -  =  ',E10.3/ 

.  'SKELETON  POISSONS  RATIO  -  =  ',F10.2/ 

•POROSITY  -  =  ',F10.2/ 

.  'COEFFICIENT  OF  PERMEABILITY  -  =  ',E10.3///) 


RW  =  0.0361 
PK  =  PK/RW 


IF(NEF.EO.l)  GO  TO  10 
READ(6,1004)  RKW,RKS 

1004  FORMAT ( 2E1 0.3) 

WRITE (7, 2004)  RKW,RKS 

2004  FORMAT('BULK  MODULUS  OF  PORE  FLUID  -  = 

.  'BULK  MODULUS  OF  SOLID  -  = 

ALPHA  =  ( POR+ ( 1 . -POR ) * ( RKW/RKS ) ) /RKW 
IF(NDF.EO.O)  GOTO  5 
BKS=CM*(1.0+V)/(3.0*(1.0-V) ) 

ALPHA=ALPHA-BKS/( RKS*RKS ) 

RSG=BKS/RKS 
5  CONTINUE 
GO  TO  20 
10  CONTINUE 

READ{6,1005)  SO,STAW 

1005  FORMAT(2F10.0) 

WRITE(7,2005)  SO,STAW 

2005  FORMAT (' INITIAL  DEGREE  OF  SATURATION  -  = 

.  'SURFACE  TENSION  OF  AIR-WATER  MIXTURE  = 

20  CONTINUE 


',E10.3/ 

',E10.3///) 


'  ,F10.2/ 

'  ,F10.2///) 


INPUT  FOR  TIME  VARIABLES 


READ(6,1006)  TEND 
1006  FORMAT(E10. 3) 

WRITE (7, 2006)  TEND 

2006  FORMAT ( 'MAXIMUM  CALCULATION  TIME  -  =  ',E10.3/) 

C 

IF(NTF.EO.l)  GO  TO  30 
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READ(6,1007)  DT 

1007  FORMAT (ElO. 3) 

WRITE{7,2007)  DT 

2007  FORMAT (  'CONSTANT  TIME  STEP - -  ',E10.3/) 

GO  TO  40 
30  CONTINUE 

READ(6,1008)  NDT 

1008  FORMAT (15) 

WRITE(7,2008)  NDT 


2008  FORMAT('NO.  OF  VARIABLE  TIME  STEPS - -  ',15  /) 

WRITE(7,2009) 

2009  FORMAT('NO.  OF  CYCLES  TIME  STEP  '  ) 


DO  35  1=1, NDT 

READ(6,1010)  NCL(I),DTT(I) 
1010  FORMAT(I5,E10.3) 

WRITE (7,2010)  NCL(I),DTT(I) 

2010  FORMAT(I5,10X,E10.3) 

35  CONTINUE 

40  CONTINUE 

WRITE(7,2011) 

2011  FORMAT (///) 

INPUT  FOR  INITIAL  STRESSES 


READ(6,1012)  SRI,STI,PI 
1012  FORMAT(3E10.3) 

WRITE( 7,2012)  SRI,STI,PI 

2012  FORMAT ( 'INITIAL  EFFECTIVE  RADIAL  STRESS  -  =  ',E10.3/ 

.  'INITIAL  EFFECTIVE  CIRCUMF.  STRESS  =  ',E10.3/ 

.  'INITIAL  PORE  FLUID  PRESSURE - : - -  ',E10.3///) 


INPUT  FOR  BOUNDARY  CONDITIONS  AT  LEFT  END 


READ( 6,1013)  JSDL, JDCL,TNFL,BPPL,SSDL 
1013  FORMAT(2I5,3E10.3) 

WRITE( 7,2013)  JSDL 

2013  FORMAT ( 'SKELETON  DISPL.  AT  LEFT  END  BOUNDARY  '/ 

.  '.EO.O  FREE  '/ 

:  '.EO.l  FIXED  OR  SPECIFIED  -  =  ',15//) 

IF(JSDL.EO.O)  WRITE(7,2014)  TNFL 

2014  FORMAT ( 'APPLIED  TOTAL  STRESS  AT  LEFT  END  =  ',E10.3//) 
IF( JSDL. EO.l)  WRITE(7,2114)  SSDL 

114  FORMAT ( 'SPEC.  SKELETON  DISP.  AT  LEFT  END  =  ',E10.3//) 
WRITE(7,2015)  JDCL 

2015  FORMAT ( 'DRAINAGE  CONDITION  AT  LEFT  END  BOUNDARY  '/ 

.  '.EO.O  PERMEABLE  '/ 

.  '.EO.l  IMPERMEABLE  -  =  ',15//) 

IF( JDCL.EO. 0)  WRITE(7,2016)  BPPL 

2016  FORMAT ( 'SPECIFIED  PORE  PRESSURE  AT  LEFT  END  =  ',E10.3//) 


INPUT  FOR  BOUNDARY  CONDITION  AT  RIGHT  END 
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READ( 6, 1017)  JSDR,JDCR,TNFR,BPPR,SSDR 
1017  FORMAT(2I5,3E10. 3) 

WRITE! 7,2017)  JSDR 

2017  FORMAT! 'SKELETON  DISPL.  AT  RIGHT  END  BOUNDARY 

.  '.EO.O  FREE 

.  '.EO.l  FIXED  OR  SPECIFIED  - 

IF! JSDR. EO.O)  WRITE!7,2018)  TNFR 

2018  FORMAT! 'APPLIED  TOTAL  STRESS  AT  RIGHT  END 
IF!JSDR.E0.1)  WRITE!7,2118)SSDR 

2118  FORMAT! 'SPEC.  SKELETON  DISP.  AT  RIGHT  END 
WRITE!7,2019)  JDCR 

2019  FORMAT! 'DRAINAGE  CONDITION  AT  RIGHT  END 

.  '.EO.O  PERMEABLE 

.  '.EO.l  IMPERMEABLE  - 

IF! JDCR. EO.O)  WRITE!7,2020)  BPPR 

2020  FORMAT! 'SPECIFIED  PORE  PRESSURE  AT  RIGHT  END 

INPUT  TIME  OF  PRINTS 
READ! 6,*)  NPRNT 

READ! 6,*)  (TPRNT(NJ) ,NJ=1,NPRNT) 

NPRNT  =  NPRNT  +  1 
TPRNT! NPRNT)  =  l.E+01 
NPRNT  =  1 


CALCULATE  INTERNAL  FLAG 

IF!NKF. EO.O. AND. NEF. EO.O. AND. NTF. EO.O)  NLNR  = 
IP!NKF.E0. 1. AND. NEF. EO.O. AND. NTF. EO.O)  NLNR  = 
IF!NKF. EO.O. AND. NEF. EO.l. AND. NTF. EO.O)  NLNR  * 
IF!NKF. EO.O. AND. NEF. EO.O. AND. NTF. EO.l)  NLNR  = 
IF!NKF. EO.l. AND. NEF. EO.l. AND. NTF. EO.O)  NLNR  = 
IF{NKF.EO. 1. AND. NEF. EO.O. AND. NTF. EO. 1)  NLNR  = 
IF!NKF. EO.O. AND. NEF. EO. 1. AND. NTF. EO. 1 )  NLNR  = 
IF!NKF. EO.l. AND. NEF. EO.l. AND. NTF. EO.l)  NLNR  = 

C 

NPK  =  0 
NPE  =  0 
NPH  =  0 
C 

IF!NLNR.E0.2.0R.NLNR.E0.3.0R.NLNR.EO. 6 )  NPK  = 
IF!NLNR.E0.1.0R.NLNR.E0.3.0R.NLNR.E0.5)  NPE  = 
IF!NLNR.E0. 1.OR.NLNR.E0.2.OR.NLNR.E0.4)  NPH  = 
RETURN 
END 
C 

SUBROUTINE  NODEG!X) 

REAL  N,L 
DIMENSION  X!l) 

COMMON/GINPT /TITLE ! 20 ) , RI ,RO,NUMEL,NKF ,NEF ,NTF 
COMMON/GPARM/MM ! 7 ) , AP , N VF , NDF 
N  =  FLOAT !NUMEL) 


V 

'/ 

M5//) 

,E10.3//) 

,E10.3//) 

V 

V 

M5//) 

=  ',E10.3///) 


I  PLANE 


non  n  o  o  o 


L  *  RO-RI 
AX  =  RI 
X(l)  =  AX 
C 

DO  100  I=1,NUMEL 

A1  =  L*(1.+(N-1. )*(4,*AP-2. )/N)/N 
B  =  (2.*L/N-2.*A1)/(N-1. ) 

AI  *  A1+(I-1. )*B 
AX  =  AX+AI 
X(I+1)  =  AX 
100  CONTINUE 
C 

RETURN 
END 

INSERT  ELEMT 

SUBROUTINE  ELEMT ( ID, I DH , IHD, IXS , IXF , IXH ,NUMEL, NUMNP , 

*  NE0,NE0H,X) 

DIMENSION  ID(NUMNP,2) ,IDH(1) ,IHD(1) ,IXS(2,1) ,IXF(2, 1) , 

*  IXH(2,1),X(1) 

COMMON/BPRES/NMOD,NOD( 2) ,SPR( 2) 

COMMON/BFORC/NMOS , NOS ( 2 ) , TNF ( 2 ) , NMOU , NOU ( 2 ) , SPU ( 2 ) 
COMMON/BCONC/JSDL, JDCL,TNFL,BPPL, JSDR, JDCR,TNFR,BPPR 

#,SSDL,SSDR 

CALL  IZERO( ID,NUMNP*2) 

CALL  IZERO(rDH,NUMNP) 

ASSIGN  SKELETON  DISPL,  BOUNDARY  CONDITIONS 

NMOS  =  0 
NMOU  =  0 

IF( JSDL.NE.O)  GO  TO  10 
NMOS  =  NMOS+1 
ID{ 1,1 )  =  0 
NOS (NMOS)  =  1 
TNF (NMOS)  =  -TNFL*X(1) 

GO  TO  20 
10  CONTINUE 
NMOU=NMOU+l 
ID(1,1)  =  0 
NOU (NMOU )=1 
SPU ( NMOU )=SSDL 
20  CONTINUE 

IF( JSDR.NE.O)  GO  TO  30 
NMOS  =  NMOS+1 
ID(NUMNP,1)  =  0 
NOS (NMOS)  =  NUMNP 
TNF(NMOS)  =  TNFR*X(NUMNP) 

GO  TO  40 
30  CONTINUE 
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NM0U»NM0U+1 
ID(NUMNF,1)  «  0 
NOU(NMOU)=NUMNP 
SPU(NMOU)=SSDR 
40  CONTINUE 

ASSIGN  DRAINAGE  BOUNDARY  CONDITIONS 
NMOD  =  0 

IF( JDCL.NE.O)  GO  TO  48 
IF(BPPL.NE.0,0)  GO  TO  42 
ID(1,2)  =  1 
IDH(l)  =  1 
GO  TO  50 
42  CONTINUE 

NMOD  =  NMOD+1 
ID(1,2)  =  0 
IDH(l)  =  0 
NOD (NMOD)  =  1 
SPR(NMOD)  =  BPPL 
GO  TO  50 
48  CONTINUE 
ID(1,2)  =  0 
IDH(l)  =  0 
50  CONTINUE 

IF(JDCR.NE.O)  GO  TO  58 
IF(BPPR.NE.O.O)  GO  TO  52 
ID(NUMNP,1)  =  1 
IDH(NUMNP)  =  1 
GO  TO  60 
52  CONTINUE 

NMOD  =  NMOD+1 
ID(NUMNP,2)  =  0 
IDH(NUMNP)  =  0 
NOD(NMOD)  =  NUMNP 
SPR(NMOD)  =  BPPR 
GO  TO  60 
58  CONTINUE 

ID(NUMNP,2)  =  0 
IDH(NUMNP)  =  0 
60  CONTINUE 

NEO  =  0 
NEOH  =  0 
DO  75  1=1, NUMNP 
DO  70  J=l,2 
IF(ID(I,J))  68,65,68 
65  NEO  =  NEO+1 
ID(I,J)  =  NEO 
IF(J.E0.2)  GO  TO  66 
GO  TO  70 
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66  NEOH  =  NEQH+l 
IDH(I)  =  NEOH 
IHD(NEOH)  =  NEO 
GO  TO  70 
68  ID(I,J)  =  -1 

IF(J.E0.2)  IDH(I)  =  -1 
70  CONTINUE 
75  CONTINUE 

DO  100  I=1,NUMEL 
IXS(1,I)  =  ID(I,1) 

IXS(2,I)  =  ID(I+1,1) 

IXF(1,I)  =  ID{I,2) 

IXF(2,I)  =  ID(I+1,2) 

IXH(1,I)  =  IDH(I) 

IXH(2,I)  =  IDH(I+1) 

100  CONTINUE 

IF(NMOS.EO.O)  GO  TO  106 
DO  104  I=l,NMOS 
II  =  NOS (I) 

104  NOS(I)  =  ID(II,1) 

106  CONTINUE 

IF(NMOD.EO.O)  GO  TO  120 
DO  110  I=l,NMOD 
II  =  NOD(I) 

110  NOD(I)  =  ID(II,2) 

120  CONTINUE 

IF(NMOU.EO.O)  GOTO  140 
DO  130  I=l,NMOU 
II=NOU(I) 

130  NOU(I)=ID(II,l) 

140  CONTINUE 

RETURN 

END 

SUBROUTINE  KSTAR ( X , GK , GKO , MB , GH , GHO , MBH , PNO , RS , I XS , I XF , I XH , 

*  S , STRES , NUMEL , NUMNP , NEQ , NEQH , NEOH  2 , SHA , PN I ) 
DIMENSION  X{NUMNP) ,GK(NE0,4) ,GKO(NEO,4) ,MB(NE0) , 

*  GH(NE0H,2) ,GHO(NEOH,2) ,MBH(NE0H2) , PNO ( NEQH ) , RS ( NEO ) , 

*  IXS(2,NUMEL) ,IXF(2,NUMEL) ,IXH(2,NUMEL) ,S(2,2,NUMEL) , 

*  STRES ( 2 , NUMEL ) , CT ( 2 , 2 ) , SHA ( NEQ ) , PP ( 2 ) , DI ( 2 ) , PN I { NEQH ) 

DIMENS ION  E(2,2),H(2,2),C(2,2),EHP(2,2), EHM ( 2 , 2 )  ,  R ( 2 ) 
COMMON/GINPT /TITLE ( 20 ) ,RI ,RO,NUMEX,NKF ,NEF ,NTF , IPLANE 
COMMON/GPARM/NLNR , NPK , NPE , NPH , NUMNX , NEX , NEQX , AP , N VF , NDF 
COMMON/MTPRT/CM , V , POR , PK , ALPHA , SO , STAW , RSG 
COMMON/TIMEV/DT,TEND,NDT,DTT( 10) ,NCL( 10) 
COMMON/INIST/SR,ST,PI 

COMMON/BPRES/NMOD,NOD( 2) ,SPR( 2) 

COMMON/BFORC/NMOS , NOS ( 2 ) ,TNF ( 2 ) , NMOU , NOU ( 2 ) , SPU ( 2 ) 
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IF(NTF.EO.l)  DT  =  DTT(l) 

CALCULATE  INITIAL  VALUE  OF  ALPHA 

IF{NEF.NE.l)  GO  TO  50 
PAO  =  14.7 
HC  =  0.02 
P  =  PAO+STAW+PI 
P2  =  P*P 

ALPHA  =  POR* (1.0-SO+HS*SO)*PAO/P2 
50  CONTINUE 

NE04  =  NE0*4 
CALL  SZERO(GK,NE04 ) 

CALL  SZERO(GH,NE0H2 ) 

CALL  IZERO(MB,NEQ) 

CALL  IZERO(MBH,NEOH2) 

CALL  SZERO(RS,NE0) 

IF(NLNR.NE.O)  CALL  SZERO(GKO,NE04 ) 

IF(NLNR.GT. 1. AND.NLNR.lt. 6)  CALL  SZERO (GHO , NE0H2 ) 

DO  200  NN=l,NOMEL 
N  -  NN 

ASSEMBLE  GLOBAL  STIFFNESS  MATRIX 


CALL  CONST(X,N,A,XL) 

C 

CALL  ELSTF(CM,V,A,XL,S(1,1,N) ) 

C 

CALL  ASMBM (GK , S ( 1 , 1 , N ) , IXS ( 1 , N  > , IXS ( 1 , N ) , NEQ , 4 ) 

IF(NPK.EO.l) 

*  CALL  ASMBM(GKO,S(l,l,N),IXS(l,N),IXS(l,N),NEO,4) 
C 

CALL  ELEMX(ALPHA,A,XL,E) 

C 

CALL  ELHMX(PK,A,XL,H) 

C 

DO  150  1=1,2 

DO  150  J=l,2 

EHP(I,J)  =  0.5*DT*H(I,J)+E(I,J) 

EHM(I,J)  =  -0.5*DT*H(I,J)+E(I,J) 

150  CONTINUE 
C 

CALL  ASMBM(GK,EHM,IXF(1,N),IXF(1,N) ,NE0,4) 

IF(NPE.EO. 1) 

*  CALL  ASMBM(GKO,E,IXF(l,N),IXF(l,N),NEQ,4) 

C 

CALL  ELCMX{A,XL,C,CT) 

IF(NDF.EO.O)  GOTO  157 

DO  155  1=1,2 
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DO  155  J*l,2 

C(I,J)=(1.0-RSG)*C(I,J) 

CT(I,J)*(1.0-RSG)*CT(I,J) 

155  CONTINUE 
157  CONTINUE 
C 

CALL  ASMBM(GK,C,IXS(l,N),IXF(i,N) ,NE0,4) 
IF(NLNR.NE,0) 

*  CALL  ASMBM(GKO,C,IXS(l,N),IXF(l,N),NEQ,4) 
IF(NLNR.NE.O) 

*  CALL  ASMBM(GKO,CT,IXF(l,N) ,IXS(1,N) ,NE0,4) 
C 

CALL  ASMBM(GK,CT,IXF( 1,N) ,IXS(1,N) ,NE0,4 ) 

C 

IF(NPH.NE.l)  GO  TO  180 
DO  160  1=1,2 
DO  160  J=l,2 

EHP(I,J)  =  0.5*DT*H(I,J) 

160  EHM(I,J)  =  -0.5*DT*H(I,J) 

CALL  ASMBM(GKO,EHM,IXF(l,N) ,IXF(1,N) ,NE0,4) 
180  CONTINUE 


ASSEMBLE  GLOBAL  GH  MATRIX 

CALL  ASMBM(GH,EHP,IXH(1,N) ,IXH(1,N) ,NEOH,2) 
IF(NLNR.LE. 1.0R.NLNR.GE.6)  GO  TO  190 
IF(NPE.EO.l) 

*  CALL  ASMBM(GHO,E,IXH(l,N),IXH{l,N),NE0H,2) 
IF(NPH.EO.l) 

*  CALL  ASMBM(GHO,EHP,IXH(l,N) ,IXH(1,N) ,NE0H,2) 

RESIDUAL  LOAD  VECTOR  DUE  TO  INITIAL  EFFECTIVE  STRESS 

190  IF(SR.EO.O.O.AND.ST.EO.O.O)  GOTO  200 
CALL  RESLS(SR,ST,A,XL/R) 

CALL  ASMBV(RS,R,IXS( 1,N) ,2) 

200  CONTINUE 


TRIANGLIZE  GLOBAL  STIFFNESS  MATRIX 

IF(NMOD.EO.O.AND.NMOU.EO.O)  GO  TO  250 
CALL  SZERO(SHA,NEO) 

IF(NMOD.EO.O)  GOTO  220 
DO  210  I=l,NMOD 
210  CALL  MODFY(NOD( I ) ,SPR( I ) ,GK,SHA,NEO,4) 
220  IF(NMOU.EO.O)  GOTO  250 
DO  230  I=l,NMOU 

230  CALL  MODFY ( NOU{ I ) ,SPU( I ) ,GK,SHA,NE0, 4 ) 


Dll 


250  CONTINUE 


C 

CALL  TRIA(NE0,4,GK,MB) 

C 

C  PROFILE  GLOBAL  GH  MATRIX 

C 

CALL  PRFIL(GH,MBH,NE0H,2) 

C 

C  ASSIGN  INITIAL  PORE  PRESSURE 

C 

DO  300  I=1,NE0H 
300  PNO(I)  =  PI 
C 
C 

C  RESIDUAL  LOAD  VECTOR  DUE  TO  INITIAL  PORE  PRESSURE 

C 

C 

IF(NDF.EQ.O.OR.PI.EO.O)  GOTO  348 
DO  346  N=1,NUMEL 

CALL  CONST(X,N,A,XL) 

CALL  ELCMX(A,XL,C,CT) 

DO  342  1=1,2 
DO  342  J=l,2 

342  C(I,J)=RSG*C(I,J) 

CALL  EXTRT(PP,IXH(l,N),PNO,2) 

DI ( 1 ) =C ( 1 , 1 ) *PP ( 1 ) +C ( 1 , 2 ) *PP ( 2 ) 
DI(2)=C(2,1)*PP(1)+C(2,2)*PP(2) 

CALL  ASMBV(RS,DI,IXS(1,N) ,2) 

346  CONTINUE 
348  CONTINUE 

IF(NVF.E0. I.AND.PI.NE.0,0)  GOTO  350 
GOTO  360 

350  CALL  COPY(PNI ,PNO,NE0H) 

360  CONTINUE 
C 

C  ASSIGN  INITIAL  EFFECTIVE  STRESSES 

C 

DO  400  I=1,NUMEL 
STRES(1,I)  =  SR 
400  STRES(2,I)  =  ST 
RETURN 
END 
C 

C  INSERT  SOLVE 
C 

SUBROUTINE  SOLVE ( IHD, X ,GK, GKO, MB , GH , GHO , MBH , PNO , PNC , DISP , 

*  U I , RS V, RS , IXS , IXF , IXH , S , STRES , NUMEL , NUMNP , NEQ , NEQH , NEQH  2 , 

*  SHA,PNI) 

INTEGER  BDUMMY 

DIMENSION  IHD(1),X(1),GK(NEO,4),GKO(NEQ,4) ,MB(1) , 

*  GH(NE0H,2) ,GHO(NEOH,2) ,MBH(1) ,PNO( 1) ,PNC( 1 ) ,DISP( 1 ) , 

*  UI(1),RSV(1),RS(1),IXS(2,1),IXF(2,1),IXH{2,1), 
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*  S(2,2,1),STRES(2,1),SHA(1),PNI(1) ,PPI(2) 

DIMENSION  DU(2),E(2,2),H(2,2),R{2),PP(2),DI(2) 
COMMON/GINPT/TITLE ( 20 ) , RI , RO, NUMEX , NKF , NEF , NTF , I  PLANE 
COMMON/GPARM/NLNR , NPK , NPE ,NPH , NUMNX , NEX , NEQX , AP , NVF , NDF 
COMMON /MTPRT/CM , V, POR , PK , ALPHA , SO , STAW , RSG 
COMMON/TIMEV/DT,TEND,NDT,DTT(10) ,NCL(10) 

COMMON /BFORC/NMOS , NOS ( 2 ) , TNF ( 2 ) , NMOU , NOU ( 2 ) , S  PU ( 2 ) 
COMMON/BPRES/NMOD , NOD ( 2 ) , SPR ( 2 ) 

COMMON/INIST/SRI , STI , PI 
C 

NFIRST  =  0 
ISTEP  =  0 
TIME  =  0.0 
0  =  0.0 
VCOR  =0.0 
C 

M  =  0 
MM  =  1 
C 

CALL  SZERO(PNC,NE0H) 

CALL  SZERO(DISP,NEO) 

CALL  SZERO(UI,NEQ) 

C 

50  CONTINUE 
VOL=0.0 
NGH  =  0 

IF{NTF.NE.l)  GO  TO  70 
M  =  M+1 

IF(MM.GT.NDT)  GO  TO  70 
IF(M.GT.NCL(MM) )  GO  TO  60 
DT  =  DTT(MM) 

GO  TO  70 
60  MM  =  MM+1 
DT  =  DTT(MM) 

M  =  1 
NGH  =  1 
70  CONTINUE 
C 

ISTEP  =  ISTEP+1 
IF(NMOU.EO.O)  GOTO  73 
IF(ISTEP.NE.2)  GOTO  73 
DO  72  1=1, NMOU 

72  SPU(I)=0.0 

73  CONTINUE 

TIME  =  TIME+DT 
CALL  SZERO(RSV,NE0) 

C 

IF(NMOS.EO.O)  GO  TO  78 
DO  74  I=l,NMOS 
II  =  NOS(I) 

74  RSV(II)  =  TNF(I) 

78  CONTINUE 
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C 

WRITE( 1,2000)  ISTEP 
2000  FORMAT (*  CYCLE  NUMBER  =  *,15) 

C 

IF(NFIRST.EO.O)  GO  TO  340 
IF(NLNR.EO.O)  GO  TO  80 
IF(NLNR.EO.3.AND.NGH.EO.0)  GOTO  80 
CALL  COPY(GK,GKO,NEO*4) 

IF(NLNR.GE. 6)  CALL  SZERO( GH ,NE0H* 2 ) 

IF(NLNR.GT. 1, AND.NLNR.lt. 6)  CALL  COPY (GH ,GHO , NEQH* 2 ) 
80  CONTINUE 

ELEMENT  PACKAGE 

DO  330  N=1,NUMEL 

CALL  EXTRT(DU,IXS( 1,N) ,UI,2) 

IF(NVF.NE.l)  GOTO  90 
CALL  EXTRT(DI,IXS(1,N) ,DISP,2) 

CALL  EXTRT(PP,IXH(1,N) ,PNO,2) 

IF(PI.EO.O.O)  GOTO  85 
CALL  EXTRT(PPI,IXH(1,N) ,PNI,2) 

PP(1)=PP(1)-PPI(1) 

PP(2)=PP(2)-PPI(2) 

5  CONTINUE 

CALL  CONST(X,N,A,XL) 

CALL  FLVOL ( DI , PP , POR , ALPHA , A , XL , VOL , I PLANE , RSG , N DF ) 
WRITE(50,*) 'N=' ,N, '  VOL*' , VOL 
WRITE ( 50 , * ) • DI ( 1 ) a ' , DI ( 1 ) , •  DI ( 2 ) = ' , DI ( 2 ) 

WRITE(50,*) 'PP(1)=' ,PP(1), •  PP{2)=' ,PP(2) 

WRITE(50,*) 'A=' ,A,  '  XL=.',XL 
WRITE(50,475) 

475  FORMAT(/) 

90  CONTINUE 

IF(NLNR.EO.O)  GO  TO  250 
IF(NLNR.E0.3.AND.NGH.E0.0)  GOTO  250 
CALL  CONST(X,N,A,XL) 

NONLINEAR  IN  SOLID  SKELETON  (NKF  =  1) 

IF(NKF.NE.l)  GO  TO  100 

UPDATE  CONSTRAINED  MODULUS (CM)  AND  POISSION  RATIO(V) 
BASED  ON  EFFECTIVE  STRESSES 

CALL  MODEL{CM,V,STRES(l,N),STRES(2,N) ) 

UPDATE  ELEMENT  TANGENT  STIFFNESS  MATRIX 

CALL  ELSTF(CM,V,A,XL,S(1,1,N) ) 

C 
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ASSEMBLE  ELEMENT  TANGENT  STIFFNESS  MATRIX 
CALL  ASMBM(GK,S(1,1,N) ,IXS(1,N) ,IXS(1,N) ,NE0,4) 
UPDATE  ELEMENT  EFFECTIVE  STRESSES 
CALL  EFFST{DU,A,XL,CM,V,STRES( 1,N) ,STRES(2,N) ) 
NONLINEAR  IN  FLUID  COMPRESSIBILITY  (NEF  =  1) 

100  IF(NEF.NE.l)  GO  TO  200 

UPDATE  ALPHA  BASED  ON  OLD  PORE  PRESSURE 
CALL  VCOMP ( PNO , IXH , SO , POR , STAW , ALPHA ) 

UPDATE  ELEMENT  E  MATRIX 
CALL  ELEMX{ALPHA,A,XL,E) 

ASSEMBLE  E  MATRIX  TO  GK 

CALL  ASMBM(GK,E,IXF( 1,N) ,IXF(1,N) ,NE0,4) 
ASSEMBLE  E  MATRIX  TO  GH 

CALL  ASMBM(GH,E,IXH( 1,N) ,IXH(1,N) ,NE0H,2) 

VARIABLE  TIME  STEPS  (NTF  =  1) 

200  IF(NTF.NE.l)  GO  TO  250 

CALCULATE  ELEMENT  H  MATRIX 

CALL  ELHMX{PK,A,XL,H) 

ASSEMBLE  -0.5*DT*H  TO  GK 

CALL  FACTR(H,4,-0.5*DT) 

DO  202  1=1,2 
DO  202  J=l,2 

202  H(I,J)  =  -.5  *  DT  *  H(I,J) 

CALL  ASMBM(GK,H,IXF(1,N) ,IXF(1,N) ,NE0r4) 

ASSEMBLE  0.5*DT*H  TO  GH 

CALL  FACTR(H,4,-1.0) 

DO  204  1=1,2 
DO  204  J=l,2 
204  H(I,J)  =  -H(I,J) 

CALL  ASMBM(GH,H,IXH(1,N) ,IXH(1,N) ,NE0H,2) 

C 
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ASSEMBLE  ELEMENT  RESIDUAL  LOAD  VECTOR  DUE  TO 
EFFECTIVE  STRESSES 

250  CONTINUE 

R(l)  =  S(1,1,N)*DU(1)+S(1,2,N)*DU(2) 

R(2)  =  S(2,1,N)*DU(1)+S(2,2,N)*DU(2) 

CALL  ASMBV(RS,R,IXS( 1,N) ,2) 

330  CONTINUE 

340  CONTINUE 


DO  350  1=1, NEO 
350  RSV(I)  =  RSV(I)-RS{I) 

C 

CALL  MLTPL(GH,PNC,PNO,MBH,NEQH,NEOH2, 2) 

CALL  ASMBV(RSV,PNC,IHD,NE0H) 

C 

IF(NFIRST.EO.O)  GO  TO  380 

CALL  OUTPT ( DISP , X , PNO , STRES , PK , NUMEL , NKF , I PLANE , 0 . TIME , DT , 
*  IXH(1,1) ,NVF,VOL,IXS,VCUR) 

380  CONTINUE 
C 

IF(NFIRST.EO.O)  GOTO  386 
IF(NLNR.EO.O)  GO  TO  386 
IF(NLNR.EQ.3.AND.NGH.E0.0)  GOTO  386 
IF(NMOD.EO.O.AND.NMOU.EO.O)  GO  TO  385 
CALL  S2ERO(SHA,NE0) 

IF(NMOD.EQ.O)  GOTO  383 
DO  382  I=l,NMOD 

382  CALL  MODFY(NOD(I),SPR(I),GK,SHA,NEO,4) 

383  IF(NMOU.EO.O)  GOTO  385 
DO  384  I=l,NMOU 

384  CALL  MODFY(NOU(I) ,SPU(I) ,GK,SHA,NE0,4) 

385  CONTINUE 
C 

CALL  TRIA(NEO,4,GK,MB) 

C 

386  CONTINUE 
C 

IF(NMOD.EO.O.AND.NMOU.EO.O)  GOTO  395 
DO  387  1=1, NEO 

387  RSV(I)  =  RSV(I)+SHA(I) 

IF(NMOD.EO.O)  GOTO  389 
DO  388  I=l,NMOD 

II  =  NOD(I) 

388  RSV(II)  =  SPR(I) 

389  CONTINUE 
IF(NMOU.EO.O)  GOTO  395 
DO  392  I=l,NMOU 

II=NOU(I) 

392  RSV(II)=SPU(I ) 
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395  CONTINUE 
C 

CALL  BACKS (NEO,GK,RSV, MB) 

C 

CALL  SZERO(PNO,NE0H ) 

CALL  SZERO(PNC,NE0H) 

KK  =  1 
C 

DO  405  1=1, NEO 
UI(I)  =  RSV(I) 

IF(IHD(KK)  .NE.I  )  GO  TO  4T30 
DISP(I)  =  UI{I) 

PNO(KK)  =  UI(I) 

KK  =  KK+1 
GO  TO  405 

400  DISP(I)  =  DISP(I)+UI(I) 

405  CONTINUE 
C 

IF(NFIRST.EO.l)  GO  TO  420 
WRITE(1,2001) 

2001  FORMAT ( '  INITIAL  PORE  PRESSURE  RESPONSE’) 
DO  410  I=1,NE0H 

410  WRITE(1,2002)  I,PNO(I) 

2002  FORMAT (1 5, ElO. 3) 

C  PAUSE  55 

WRITEd,*)  '  PAUSE  55' 

READ( 1,2022)  BDUMMY 
2022  FORMAT (A2) 

IF(BDUMMy.E0.84)  STOP  84 


420  NFIRST  =  1 

IF(TIME.LT.TEND)  GO  TO  50 

RETURN 

END 

SUBROUTINE  CONST ( X , N , A , L) 

REAL  L,LL,X(1) 

COMMON/CNSTS/BO , B ( 1 2 ) 

A  =  X(N) 

L  =  X(N+1)-X(N) 

AA  =  A*A 
LL  =  L*L 

BO  =  0.5*L*(2.*A+L) 

B(l)  =  L 

B(2)  =  LOG(l.+L/A) 

B(3)  =  B(1)-A*B(2) 

B(4)  =  L/(A*(A+L)) 

B(5)  =  B(2)-A*B(4) 

B(6)  =  B(1)-2.*A*B(5)-AA*B(4) 

B(7)  =  0.5*LL 

B(8)  =  B0-2.*A*B(3)-AA*B(2) 

B(9)  =  B0-3.*A*B(6)-3.*AA*B(5)-AA*A*B(4) 
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B(10)  =*  LL*L/3. 

RINSD  =  L**3+3.*AA*L+3.*A*LL 

B(ll)  =  (RINSD)/3.-3.*A*B(8)-3.*AA*B(3)-AA*A*B(2) 

B(12)  =  0.25*(L**4) 

RETURN 
END 

ELEMENT  STIFFNESS  MATRIX 

SUBROUTINE  ELSTF ( CM , V, A, L,S ) 

REAL  L,LL,S{2,2) 

COMMON/CNSTS/BO,Bl,B2,B3,B4,B5,B6,B7,B8,B9,BlO,Bll,Bl2 
D  =  CM*V/(1.-V) 

G  =  0.5*CM*(1.-2.*V)/(1.-V) 

DP  =  D+2.*G 
LL  =  L*L 

511  =  A*DP*B1/LL-2.*A*D*B2/L+(2.*A*D/LL-2.*D/L)*B3+A*DP*B4 
+(DP-2. *A*DP/L) *B5+(A*DP/LL-2. *DP/L ) *B6+DP*B7/LL 
+2. *D*B8/LL+DP*B9/LL 

512  =  -A*DP*Bl/LL+A*D*B2/L+(D/L-2. *A*D/LL)*B3+A*DP*B5/L 
+ ( DP/L-A*DP/LL ) *B6-DP*B7/LL-2 . *D*B8/LL-DP*B9/LL 

S22  =  A*DP*B1/LL+2. *A*D*B3/LL+A*DP*B6/LL+DP*B7/LL 
+2 . *D*B8/LL+DP*B9/LL 

S(l,l)  =  Sll 
S(l,2)  =  S12 
S(2,l)  =  S12 
S(2,2)  =  S22 
RETURN 
END 

ELEMENT  C  MATRIX 

SUBROUTINE  ELCMX( A , L , C ,CT ) 

REAL  L,LL,C(2,2) ,CT(2,2) 

COMMON/CNSTS/BO,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12 
LL  =  L*L 
C 

Cll  =  -A*B1/L+(A/LL-1./L)*B7+A*B2+{1.-2.*A/L)*B3 
.  +(A/LL-2./L)*B8+B10/LL+B11/LL 

C 

C12  =  -A/LL*B7+A*B3/L+(1./L-A/LL)*B8-B10/LL-B11/LL 
C 

C21  =  A*B1/L+(1./L-A/LL)*B7+A*B3/L+(1./L-A/LL)*B8 
.  -BlO/LL-Bll/LL 

C 

C22  =  A*B7/LL+A*B8/LL+B10/LL+B11/LL 
C 

C(l,l)  =  Cll 


D18 


ono  noon  nno 


C(l,2)  =  C12 
C(2,l)  =  C21 
C(2,2)  =  C22 
C 

CT(1,1)  =  C(l,l) 

CT(1,2)  *  C(2,l) 

CT(2,1)  =  C(l,2) 

CT(2,2)  =  C(2,2) 

C 

RETURN 
END 

ELEMENT  H  MATRIX 

SUBROUTINE  ELHMX ( K , A, L, H ) 

REAL  L,LL,K,H(2,2) 

COMMON/CNSTS/BO,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12 
LL  =  L*L 

Hll  =  A*K*B1/LL+K*B7/LL 

H12  =  -Hll 

H22  =  Hll 

H(l,l)  =  Hll 
H(l,2)  =  H12 
H(2,l)  =  H12 
H(2,2)  =  H22 
RETURN 
END 

ELEMENT  E  MATRIX 

SUBROUTINE  ELEMX( ALPHA, A,L,E ) 

REAL  L,LL,E(2,2) 

COMMON/CNSTS/BO ,Bl,B2,B3,B4,B5,B6,B7,B8,B9,Bl0,Bll,Bl2 
LL  =  L*L 
C 

Ell  =  A*Bl+(l,-2.*A/L)*B7+(A/LL-2./L)*BlO+Bl2/LL 
C 

E12  =  A*B7/L+( 1./L-A/LL)*B10-B12/LL 
C 

-  E22  =  A*B10/LL+B12/LL 
C 

E(l,l)  =  -E11*ALPHA 
E(l,2)  =  -E12*ALPHA 
E(2,l)  =  -E12*ALPHA 
E(2,2)  =  -E22*ALPHA 
RETURN 
END 
C 
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ELEMENT  RESIDUAL  EFFECTIVE  STRESS  VECTOR 

SUBROUTINE  RESLS ( SR , ST , A , L, R ) 

REAL  L,R(2) 

COMMON/CNSTS/BO,Bl ,B2,B3,B4 ,B5, B6,B7,B8,B9,B10,B11,B12 

Rl  *  SR*(-A*Bl/L-B7/L)+ST*(A*B2+( l.-A/L)*B3-B8/L) 

R2  »  SR*(A*B1/L+B7/L)+ST*(A*83/L+B8/L) 

R(l)  =  Rl 
R(2)  =  R2 
RETURN 
END 

INSERT  ALIB 

SUBROUTINE  SZERO(A,N) 

DIMENSION  A(N) 

DO  100  1=1, N 
100  A(I)  =  0.0 
RETURN 
END 

SUBROUTINE  IZERO(M,N) 

DIMENSION  M(N) 

DO  100  1=1, N 
100  M(I)  =  0 
RETURN 
END 

SUBROUTINE  COPY(A,B,N) 

DIMENSION  A(1),B(1) 

DO  100  1  =  1, N 
100  A(I)  =  B(I) 

RETURN 

END 

SUBROUTINE  FACTR ( A , N , FAC ) 

DIMENSION  A(N) 

DO  100  1=1, N 
100  A(I)  =  FAC*A(I) 

RETURN 

END 

SUBROUTINE  EXTRT( DU,NA,GU,N) 

DIMENSION  DU( 1 ) ,NA( 1 ) ,GU( 1 ) 

CALL  SZERO(DU,N) 

DO  100  1=1, N 
J  =  NA(I) 

IF(J.LE.O)  GO  TO  100 
DU(I)  =  GU{J) 

100  CONTINUE 
RETURN 
END 

SUBROUTINE  VCOMP ( PNO, IXH ,SO , FOR, STAW , ALPHA ) 
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DIMENSION  PNO(l) ,IXH(1) 

PAO  *  14.7 
HC  =  0.02 
P  =  0.0 
DO  100  1=1,2 
N  =  IXH(I) 

IF(N.LE.O)  GO  TO  100 
P  =  P-PNO(N) 

100  CONTINUE 

P  =  PAO+STAW+P/2. 

PP  =  P*P 

CW  =  (1.0-SO+HC*SO) *PAO/PP 
ALPHA  =  POR*CW 

RETURN 

END 

SUBROUTINE  MODEL ( CM , V, SR,ST ) 


RETURN 

END 

SUBROUTINE  EFFST ( DU , A , L , CM , V,SR,ST ) 

REAL  DU ( 2 ) , L 
U1  *  DU(1) 

U2  =  DU(2) 

AL  =  2.*A+L 
D  =  CM*V/(1.-V) 

G  =  0.5*CM*(1.-2.*V)/(1.-V) 

DP  =  D^2.*G 

SR  =  (D/AL-DP/L)*U1+(D/AL+DP/L)*U2 
ST  =  (DP/AL-D/L)*U1+(DP/AL+D/L)*U2 
C 

RETURN 

END 

C 

SUBROUTINE  OUTPT ( DISP , X , PNC ,STRES , K , NUMEL , NKF , iPLANE , 0 » 
*  TIME, DT,LF,NVF, VOL, IXS,VCUR) 

REAL  DISP( 1) ,X(1) ,PNC(1) ,STRES{ 2,1) ,L,K 
DIMENSION  IXS{2,1) 

COMMON/BCONC/JSDL, JDCL,TNFL,BPPL, JSDR, JDCR,TNFR,BPPR 
#,SSDL,SSDR 

COMMON/PRINT/TPRNT(20) ,NPRNT 
C 

PIE  =  3.1415926 
L  =  X(2)-X(l) 

IF(LF.GT.O)  PDF  =  PNC(2)-PNC(1) 

IF(LF.LE.O)  PDF  =  PNC { 1 ) 
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IF(JDCL.EO.l)  PDF  =  0.0 
VR  «  -K*PDF/L 
VOLD»VCUR 
VCUR*VR 

VAVG=0 . 5* ( VOLD+VCUR ) 

435  IF(IPLANE.EO.O)  FAC  *  2. *PIE* 1 . 0*X( 1 ) 

IF(IPLANE.EO.l)  FAC  *  1.0 
IF(NVF.EO.O)  0  *  0+FAC*VR*DT 
IF(NVF.EO.l)  0  =  VOL 
C 

DPL=DISP( 1 ) 

II=IXS(2,NUMEL) 

DPR»DISP(II) 

C 

TIMED  =  TIME-DT 

OUTPUT  TIME (SEC)  AND  FLOW (GALS  PER  UNIT  FOOT) 

WRITE(9,*)  TIMEO/86400. ,0/1728.  *  7.479  *  12. 

WRITE(5,*)  TIMEO/86400. ,VAVG  *2. 0*PIE*X( 1 ) *12*60/231 
NFR=NFF+1 

IF  TIMED  LESS  THAN  PRINT  TIME  SKIP  OUTPUT  SECTION 

IF(TIMEO.LT.TPRNT(NPRNT) )  GOTO  300 
NPRNT  *  NPRNT  +  1 

WRITE(7,2001)  TIMEO 
2001  FORMAT('TIME  =  ',E10.3//) 


WRITE(7,2002)  DPL,DPR,VR,Q 

2002  FORMAT ( 'DISPLACEMENT  AT  TUNNEL  SURFACE  -  =  ',E10.3/ 

.  'DISPLACEMENT  AT  REMOTE  BNDRY - -  ',E10.3/ 

,  'FLOW  VELOCITY  AT  TUNNEL  SURFACE  -  =  ',E10.3/ 

.  'ACCUMULATED  VOLUME  OP  WATER  '  / 

.  'PER  UNIT  LENGTH  OF  TUNNEL - -  ',E10.3//) 


IF(NKF.EO.O)  WRITE(7,2003) 

2003  FORMAT ( 'DISTANCE  PORE  PRESS.’) 

C  IF(NKF.EO.l)  WRITE(7,2004) 

2004  FORMAT ( 'DISTANCE  PORE  PRESS,  EFF.  RAD.  S.  EFF.  CIR.  S.') 
DO  200  I=1,NUMEL 

IF(LF.GT.O)  GO  TO  100 
P  =  PNC(I) 

R  =  X(I+1) 

GO  TO  150 

100  P  =  0.5*(PNC(I )+PNC(I+l) ) 

R  =  0.5*(X(I)+X(I+1) ) 

150  CONTINUE 

IF(IPLANE.EO.l)  R  =  R-X(l) 

C  IF(NKF.EO.O)  WRITE(7,2005)  R,P 
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IF(NKF.EQ.O)  WRITE(8,*)  R/12. ,-P 

2005  FORMAT (F10.0,5X,F10.0) 

C  IF(NKF.EO.l)  WRITE(7,2006)  R, P, STRES ( 1 , I ) , STRES ( 2 , I ) 

2006  FORMAT(4(F10.0,5X) ) 

200  CONTINUE 

C  WRITE(7,2007) 

2007  FORMAT (///) 

WRITE(8,*)  •  0.0  123456789.  TIME  =',TIMEO 

300  CONTINUE 
RETURN 
END 

INSERT  MLIB 


MATRIX  ASSEMBLE 

SUBROUTINE  ASMBM ( ST , EST , NDR , NDC , NEQ , MBAND ) 
DIMENSION  ST (NEQ, MBAND) ,EST(2,2) ,NDR( 2 ) ,NDC ( 2 ) 
DO  200  K=l,2 
IR  =  NDR(K) 

IF(IR.LE.O)  GO  TO  200 
DO  100  L=l,2 
IC  =  NDC(L) 

IF(IC.LE.O)  GO  TO  100 
IC  =  IC-IR+1 
IF(IC.LE.O)  GO  TO  100 
ST(IR,IC)  »  ST(IR,IC)+EST(K,L) 

100  CONTINUE 
200  CONTINUE 
RETURN 
END 

VECTOR  ASSEMBLE 

SUBROUTINE  ASMBV( VC , EVC ,NDR ,NA) 

DIMENSION  VC(1),EVC(1),NDR(1) 

DO  300  K=1,NA 
IR  =  NDR(K) 

IF(IR.LE.O)  GO  TO  300 
VC(IR)  =  VC(IR)+EVC(K) 

300  CONTINUE 
RETURN 
END 

TRIANGLIZE 

SUBROUTINE  TRIA( NEQ, M , A,MB ) 

DIMENSION  A(NE0,1),MB(NE0) 

NE  =  NEO-1 
MN  =  M-1 
MM  =  M 
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MK  =  MEO-MN 
DO  3L0  N»1,NE 
NT  =  N-MK 

IF(NT.GT.O)  MM  =  MM-1 
MB(N:  =  0 

IF(A(N,1) .EQ.O.O)  GO  TO  300 
L  =  N 
IH  *  MM 
JB  =  0 
IB  =  0 

DO  200  1=2, IH 
L  =  L+1 
J  =  L 
IB  =  IB+1 
AI  =  A(N,I) 

C  =  AI/A(N,1) 

IF(C.EO.O.O)  GO  TO  200 
JC  =  1 

DO  100  K=I,IH 

A(J,.;C)  =  A(  J,JC)-C*A(N,K) 

100  JC  =  JC+1 
A(N,:)  =  C 
JB  =  IB 
200  CONTINUE 
MB ( N :  =  JB 
300  CONTINUE 

MB(NKO)  =  0 

RETUi^N 

END 

BACK.-UBSTITUTE 

SUBR(  UTINE  BACKS(NN,A,B,MB) 

DIME  SION  A(NN,1) ,B(NN) ,MB(NN) 

N  = 

270  N  =  ■  +1 
C  =  r.(N) 

IF(A. N,l) .NE.0.0)  B(N)  =  B(N)/A(N,1) 
IF(N.EO.NN)  GO  TO  300 
IL  =  N+1 
IH  =  N+MB(N) 

M  =  : 

DO  2>5  I=IL,IH 
M  =  M  +  1 

285  B(I)  =  B(I)-A(N,M)*C 
GO  T'  270 
300  IL  =  N 
N  =  ;m-1 

IF(N.EO.O)  RETURN 
IH  =  N+MB(N) 

M  =  1 
C  =  3(N) 
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DO  400  I*IL,IH 
M  *  M+1 

400  C  =  C-A(N,M)*B(I) 

B(N)  =  C 
GO  TO  300 
END 

MATRIX  MULTIPLICATION 

SUBROUTINE  MLTPL( A , B , BO, MB ,NE0, NE02 , MBAND ) 
DIMENSION  A(NEQ,MBAND) ,B (NEQ) ,BO( NEQ ) ,MB (NEQ2 ) 
DO  300  N=1,NEQ 
BB  =  A(N,l)*BO(N) 

L  =  N 
NI  =  MB(N) 

IF(NI.GT,1)  GO  TO  50 
GO  TO  120 
50  DO  100  M=2,NI 
L  =  L+1 

100  BB  =  BB+A(N,M)*BO(L) 

120  L  =  N 

NNEO  =  N+NEQ 
NJ  =  MB (NNEO) 

IF(NJ.GT.l)  GO  TO  150 
GO  TO  250 
150  DO  200  M=2,NJ 
L  =  L-1 
NMl  =  N-M+1 

200  BB  =  BB+A(NMl,M)*BO(L) 

250  B(N)  =  B(N)+BB 
300  CONTINUE 
RETURN 
END 

MODIFY  GLOBAL  STIFFNESS  MATRIX 

SUBROUTINE  MODF Y ( N , X , A , B , NEQ , MBAND ) 

DIMENSION  A(NE0,1),B(1) 

DO  250  M=2, MBAND 


K  = 

N-M+1 

IF(K 

.LE 

.0) 

GO 

TO 

235 

B(K) 

s 

B(K) 

-A( 

K,M 

) 

*X 

A(K, 

M) 

=  0. 

0 

235 

L  = 

N+M 

-1 

IF (NEQ- 

L)  245, 

240 

9 

240 

240 

B(L) 

S 

B(L) 

-A( 

N,M 

) 

*X 

245 

A(N, 

M) 

=  0. 

0 

250 

CONTINUE 

A(N, 

1) 

=  1. 

0 

RETURN 

END 
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PROFILE  MATRIX 

SUBROUTINE  PRFIL( A , MB , NEQ ,MBAND ) 
DIMENSION  A(NE0,1),MB(1) 

DO  300  N=1,NE0 
NI  =  1 
NJ  =  1 

DO  100  M=2,MBAND 
IF(A(N,M) .NE.0.0)  NI  =  M 
IF(N-M)  100,80,80 
80  NMl  =  N-M+1 

IF(A(NM1,M) .NE.0.0)  NJ  =  M 
100  CONTINUE 
MB(N)  =  NI 
NNEQ  =  N+NEQ 
MB(NNEO)  =  NJ 
300  CONTINUE 
RETURN 
END 


SUBROUTINE  FL VOL ( DU , PP , POR , ALPHA , A , XL , VOL , I  PLANE , RSG , N DF ) 
DIMENSION  DU(2),PP(2) 

COMMON/CNSTS/BO,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12 

AL1=1.0/XL 

AL2=1.0/(2.0*A+XL) 

EV=DU ( 1 ) * ( AL2-AL1 ) +DU ( 2 ) * ( AL2+AL1 ) 

IF(IPLANE.EO.O)  VEL  =  ( A*B1+B7 ) *2 . 0* 3 . 14 1 5926 
IF(IPLANE.EO.l)  VEL  =1,0*XL 
P  =  0.5*(PP(1)+PP(2) ) 

IF(NDF.EO.O)  VOL  =  VOL- ( EV-ALPHA*P ) *VEL 

IF(NDF.EO.l)  VOL  =  VOL- ( EV* { 1 . O-RSG ) -ALPHA*? )* VEL 

WRITE(50,*) 'VOL  IN  FLVOL= ' , VOL, '  EV=',EV,’  P= ' , P 

WRITE(50,*) 'VEL=' ,VEL, '  A=',A,’  XL=’,XL,'  Bl=',Bl,'  B7=',B7 

RETURN 

END 


D26 


REPORT  DOCUMENTATION  PAGE 

Form  Approvod 

OMB  No.  0704-0188 

Mik  ftpertiM  burPtfi  fpr  tin*  cetttcttOA  oi  information  n  ntunatM  to  \  hour  ptt  roiponM.  mcMiM  tha  urn*  for 

oatfMrMM  andmaintamma  tha  dau  naadad.  and  complattng  and  raviawinQ  tha  <oda<Uon  of  inf  ormatton.  Sand  commant*  raoardiM  thi*  attiinata  or  y*y  oovar  mpm  oftw 

coSSSon  of  mSrmaSoiOnelyding  tuMtitlona  for  radueing  thn  burdan.  lo  Wa$h«ngu»n  madouartara  Sarvicai.  Oiractorata  w  mformaiw  arid  JJapo^  12IS  Jaffanon 

OavH  tlljfwaif,  Suita  1204,  ArNn9lon7>^22202-4302.  and  10  tlia  Offka  of  Managamant  and  Oudgat.  faparurorfc  Raduction  arofact  (0704-0tM),  Washington,  DC  20503- 

1.  AGEN^  USE  ONLY  (L*»vt  blank)  2.  REPORT  DATE  I  3.  REPORT  TYPE  ANO  OATES  COVERED 

September  1994  |  Final  report 

4.  TIUE  ANO  SUBTITLE 

Task  IV:  Groundshock-Induced  Hydrogeologic  Response;  Volume  II: 
Hydrologic  Response  of  Deep  Based  Systems  to  Blast  Loading 

S.  FUNDING  NUMBERS 

DACW45-84-C-0128 

C.  AUTHORIS) 

Scott  £.  Blouin,  Kwang  J.  Kim 

7.  PERFORMING  ORGANIZATION  NAME(S)  ANO  AOORESS(ES) 

Applied  Research  Associates,  Inc. 

New  England  Division 

South  Royaltion,  VT  05068 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

Miscellaneous  Paper 

GL-94-49 

9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

See  reverse. 

10.  SPONSORING /MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

The  main  text  was  published  under  separate  cover.  Available  from  the  National  Technical  Information 

Service,  S28S  Port  Royal  Road,  SpringHeld,  VA  22161. 

12a.  DISTRIBUTION /availability  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  worth) 

Meclumisms  controlling  the  hydrologic  response  of  a  deep  based  system  to  nuclear  attack  are  identified 
through  an  analysis  of  data  from  underground  tunnel  complexes  at  the  Nevada  Test  Site.  Two  mechanisms  are 
identified,  inflow  of  fracture  water  mobilized  by  the  dynamic  loadings  and  inflow  of  pore  water  due  to 
explosion  generated  residual  stresses  and  pore  pressures.  Both  mechanisms  are  serious  potential  threats  to 
deep  based  systems  located  in  geologies  similar  to  Generic  Mountain  C.  The  magnitude  of  the  fracture  water 
thrrat  is  impossible  to  assess  due  to  nonavailability  of  sufficient  geologic  and  hydrologic  data  at  Mountain  C. 

The  magnitude  of  the  pore  water  threat  at  site  C  is  estimated  based  on  flow  calculations  using  the  one 
dimensional  axisymmetric  code  CONSL.  Calculations  indicate  that  pore  water  flow  at  the  Generic  Mountain 

C  site  will  be  orders  of  magnitude  greater  than  similar  flow  at  the  Nevada  Test  Site  because  of  the  much  more 
permeable  rock  at  Mountain  C.  The  need  for  more  refrned  calculations  of  the  Mountain  C  response  is 
emphasized. 

14.  SUBJEa  TERMS 

See  reverse. 

15.  NUMBER  OF  PAGES 

219 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 

IB.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  UMITATION  OF  ABSTRACT 

NSN  754(M)t-280-5S00 


Standard  Form  298  (Rev.  2-89) 

Prncribcd  by  ANSI  Std  ZJ9-1S 

ne-102 


9.  (Continued). 


U.S.  Anny  Engineer  Waterways  Experiment  Station 
3909  Halls  Ferry  Road 
Vicksburg.  MS  39180*6199 

U.S.  Air  Force  Ballistic  Missile  Office 
Norton  Air  Force  Base 


14.  (Continued). 

Consolidation 
Deep  basing 
Explosive  loadings 


Hydrologic  response 
Pore  pressure  dissipation 
Rock  mechanics 


